C . THE N U M E R I C A L S O L U T I O N OF A S Y S T E M OF D I F F U S I O N / C O N V E C T I O N EQUATIONS DESCRIBING COASTAL CIRCULATION by THOMAS 0. B.Sc, The U n i v e r s i t y A THESIS THE SUBMITTED Of NICOL British Columbia, IN PARTIAL FULFILMENT REQUIREMENTS MASTER FOR OF THE DEGREE 1976 OF OF SCIENCE in THE Institute We FACULTY Of A p p l i e d accept to THE this OF GRADUATE Mathematics thesis the required UNIVERSITY OF October © Thomas STUDIES 0. as And Statistics conforming standard BRITISH COLUMBIA 1983 Nicol," 1983 l In presenting requirements for Columbia, I available for permission this in shall reference and study. I extensive c o p y i n g of representatives. allowed without D e p a r t m e n t of this thesis my w r i t t e n Institute The U n i v e r s i t y of B r i t i s h 2075 Wesbrook P l a c e V a n c o u v e r , Canada V6T 1W5 Date: October It 14, 1983 for Head of is this of U n i v e r s i t y of Library or of the the may be g r a n t e d by t h e publication fulfilment that purposes her partial an a d v a n c e d d e g r e e a t agree for thesis make further thesis financial gain British it freely agree for that that scholarly my D e p a r t m e n t o r understood the by his copying shall not permission. Of A p p l i e d M a t h e m a t i c s And Columbia Statistics or be i i Abstract A mathematical examined. The derivation of equations partial and consist model of differential the problem A difference linear a r e examined. and i s d i s c u s s e d and a The scheme, is method, the presented for e q u a t i o n s , and the their The i n i t i a l convergence Methods f o r d e a l i n g w i t h effects on accuracy s p e c i a l problems, and methods A finite solve them. more efficient scheme f o r o u r e q u a t i o n s , b a s e d on t h e is these methods is developed describe difference formula, A computer program t o s o l v e t h e model e q u a t i o n s , and r e s u l t s for test with f o r a d d i t i o n s t o t h e program t o model t h e recommendations c u r r e n t s more a c c u r a t e l y . data we Crank-Nicolson i n t r o d u c e d , and i t s advantages d i s c u s s e d . incorporating and and boundary c o n d i t i o n s of t h e model e q u a t i o n s p r e s e n t to model nonlinear, parabolic numerical o f t h e method d i s c u s s e d . terms, ocean c u r r e n t s i s presented. p a i r of coupled, equations. finite nonlinear stability, for equations a the corresponding stability the motivation the Crank-Nicolson solving model d e s c r i b i n g n e a r s h o r e presented. We conclude iii T a b l e of C o n t e n t s Abstract L i s t of F i g u r e s Acknowledgements I. INTRODUCTION 11 . CHAPTER 1 - 2.1 D e r i v a t i o n Of The M o d e l 2.2 The Domain Of The M o d e l 2.3 B o u n d a r y C o n d i t i o n s 2.4 I n i t i a l C o n d i t i o n s 2.5 How The M o d e l W i l l Be U s e d III. CHAPTER 2 3.1 N o t a t i o n 3.2 F i n i t e D i f f e r e n c e s 3.3 C o n s i s t e n c y , S t a b i l i t y And C o n v e r g e n c e 3.4 V a r i a b l e C o e f f i c i e n t s 3.5 S t a b i l i t y Of B o u n d a r y C o n d i t i o n s 3.6 N o n l i n e a r Terms 3.7 S o l v i n g The E q u a t i o n s IV. CHAPTER 3 4.1 I n t r o d u c t i o n To ADI 4.2 A c c u r a c y And S t a b i l i t y 4.3 C o m p u t a t i o n a l D e t a i l s And Sample P r o b l e m s 4.4 A p p l y i n g ADI To The M o d e l E q u a t i o n s V. CHAPTER 4 5.1 P r e l i m i n a r y T e s t s 5.2 I n i t i a l C o n d i t i o n s A t The L a n d B o u n d a r y 5.3 F i n d i n g H On The B o u n d a r i e s 5.4 B o u n d a r y C o n d i t i o n F o r S 5.5 T e s t s Of The M o d e l E q u a t i o n s 5.6 I n t e r p r e t a t i o n Of R e s u l t s 5.7 R e c o m m e n d a t i o n s F o r P r o c e e d i n g ii iv v 1 4 4 8 8 9 10 13 13 14 19 24 25 27 29 32 33 35 37 42 44 44 46 49 52 53 54 57 BIBLIOGRAPHY 78 iv List 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. 29. 30. 31. 32. 33. 34. 35. of F i g u r e s Domain o f t h e C o a s t a l C i r c u l a t i o n Domain o f t h e N u m e r i c a l M o d e l Test 1 : S a t time s t e p 2 Test 1 : S a t time s t e p 5 T e s t 1 : S a t t i m e s t e p 10 T e s t 1 : S a t t i m e s t e p 15 T e s t 1 : S a t t i m e s t e p 20 T e s t 1 : S a t t i m e s t e p 30 T e s t 1 : S a t t i m e s t e p 40 T e s t 1: S a t t i m e s t e p 50 Test 1 : H a t time s t e p 2 Test 1 : H at time s t e p 5 T e s t 1 : H a t t i m e s t e p 10 T e s t 1 : H a t t i m e s t e p 15 T e s t 1: H a t t i m e s t e p 20 T e s t 1 : H a t t i m e s t e p 30 T e s t 1 : H a t t i m e s t e p 40 T e s t 1 : H a t t i m e s t e p 50 T e s t 1: H S a t t i m e s t e p 2 T e s t 1: H S a t t i m e s t e p 5 T e s t 1: H S a t t i m e s t e p 10 T e s t 1: H S a t t i m e s t e p 15 T e s t 1: H S a t t i m e s t e p 20 T e s t 1: H S a t t i m e s t e p 30 T e s t 1: H S a t t i m e s t e p 40 T e s t 1: H S a t t i m e s t e p 50 T e s t 2: S a t t i m e s t e p 20 T e s t 2: S a t t i m e s t e p 30 T e s t 2: S a t t i m e s t e p 40 T e s t 2: H a t t i m e s t e p 20 T e s t 2: H a t t i m e s t e p 30 T e s t 2: H a t t i m e s t e p 40 T e s t 2: H S a t t i m e s t e p 20 T e s t 2: H S a t t i m e s t e p 30 T e s t 2: H S a t t i m e s t e p 40 2 2 2 2 2 2 2 2 2 2 2 Model 11 12 60 61 61 62 62 63 63 64 64 65 65 66 66 67 67 68 68 69 69 70 70 71 71 72 72 73 73 74 75 75 76 77 77 V Acknowledgement The c o m p l e t i o n given over suggestions Emery; of many this thesis hours, from P a u l and from a l l , has from Leblond, a great depended Jim Varah, and t h e deal of on t h e the original patience. good advice, multitude idea of of Bill 1 I. S t u d i e s of ocean that fresh water circulation consider closely the of coastal regions i s an i m p o r t a n t f a c t o r of l a r g e r i v e r s run-off Department of at o f f the n o r t h e r n Washington State e f f e c t s on a l o n g s h o r e aspect coast of in studies ocean, the r e g i o n h a s n o t been P a u l L e B l o n d and Dr. Oceanography shown in large-scale coastal over a l a r g e c o a s t a l Dr. have While these on t h e William Emery of the U n i v e r s i t y of B r i t i s h have u n d e r t a k e n a s t u d y of r i v e r distribution One in (see [6] f o r s e v e r a l examples). examined. Columbia currents run-off the e f f e c t influence INTRODUCTION run-off and salinity of s o u t h e r n B r i t i s h Columbia and an effort to understand their currents. the study has been t h e d e v e l o p m e n t of a m a t h e m a t i c a l model t o d e s c r i b e t h o s e c o a s t a l c u r r e n t s c r e a t e d by river run-off. describe salinity, The m o d e l c o n s i s t s o f a s e t nearshore hydrodynamics in of equations terms of changes a n d t h u s c a n be u s e d t o m e a s u r e t h e i m p a c t o f run-off. The main purpose of this thesis that in coastal i s t o determine a n u m e r i c a l method t o s o l v e t h e e q u a t i o n s d e f i n i n g t h e model, and to d e v e l o p a computer program i m p l e m e n t i n g t h e method. As w e l l as b e i n g of i n t e r e s t from a p h y s i c a l p o i n t of v i e w , i t i s hoped t h a t t h e r e s u l t s be used to particular, gain an from t h e understanding of o f f s h o r e i t h a s been o b s e r v e d t h a t percentage of F r a s e r R i v e r salmon in oceanographic most model may fisheries. years a In large m i g r a t e v i a t h e S t r a i t of Juan de Fuca, w h i l e the remainder return through Johnstone the n o r t h e r n end of Vancouver Island. S t r a i t at I n some y e a r s , h o w e v e r , a 2 m a j o r i t y return through Johnstone S t r a i t . are not known, but i t is likely c i r c u l a t i o n plays a large part understanding this of these m e c h a n i s m s may route is chosen. An l e a d t o an e x p l a n a t i o n f o r One c o n t a i n s a d i s c u s s i o n o f t h e d e r i v a t i o n o f t h e variables initial and parameters In C h a p t e r Two, difference accuracy follow, and equations, specified stability equations discussed. The examined. relate a linearized boundary are we by the how the examine finite consistency, form of their the effects terms i n t h e model the on equations presented. systems Crank-Nicolson More e f f i c i e n t of linear scheme, a n d d i s c u s s difference the a l t e r n a t i n g d i r e c t i o n They a r e d e s c r i b e d equation. of c o n d i t i o n s a r e then d i f f e r e n c e s and The n o n l i n e a r Three, produced b o u n d a r y and t o a p p l i c a t i o n of the and methods o f d e a l i n g w i t h them Chapter schemes implicit, for or ADI, a n d t h e c o n c e p t s o f C h a p t e r Two a r e t h e P e a c e m a n - R a c h f o r d ADI scheme, for Methods of t r e a t i n g t h e a d d i t i o n a l problems by t h e m o d e l e q u a t i o n s , derivatives, description The c o n c e p t s o f they of e x t e n d e d t o one o f t h e s e , linear the of the chapter. differencing methods t o s o l v e them. methods. as i n terms of f i n i t e equations and Definitions t h e n o t a t i o n t o be u s e d t o d e s c r i b e stability, are are discussed In A schemes i s i n t r o d u c e d . Crank-Nicolson presented t o be s o l v e d . conditions are s p e c i f i e d . m o d e l i s t o be u s e d c o n c l u d e s a that the l a r g e - s c a l e c o a s t a l i n which problem, l e a d i n g to the equations the this behavior. Chapter the The r e a s o n s f o r variable including coefficients and lower nonlinear order space terms, are 3 e x a m i n e d , and tested by solutions. Results a c c u r a c y and stability discretize the solving from full of sample computer the problems runs are shown solutions discussed. e q u a t i o n s i n t e r m s of with the known and the Finally, we Peaceman-Rachford scheme. I n C h a p t e r F o u r , t h e ADI equations with simplified method i s applied of the initial conditions t h i s problem i s d e a l t with. thickness to ADI layer method conditions full model for c o n c l u d e the explained. calculating i s presented. equations interpretations. The i s e x a m i n e d and h a n d l e them a r e the b o u n d a r i e s to ensure the scheme s a t i s f i e s b a s i c p h y s i c a l p r o p e r t i e s . nature to i s discussed, difference discontinuous and we show how for the l a n d boundary c o n d i t i o n the p r o b l e m s and Next the the salinity given, methods we formula required Sample r u n s and are The model defect their with by r e s u l t s for brief the boundary the physical Recommendations f o r p r o c e e d i n g w i t h the chapter. use model 4 II. 2.1 D e r i v a t i o n Of The The CHAPTER 1 Model p r i m a r y c o n c e r n of problems presented circulation. The by this the t h e s i s i s the mathematical d e r i v a t i o n that follows justification, from for a s s u m p t i o n s and a c c o m o d a t i o n s made final equations. i m p l i c a t i o n s of A model an of this and five days, by I f we and 1, 0 the force. The We U 0 includes viewpoint, the d i s c u s s i o n of the Our derivation consider accelerations is s p a c e s c a l e , Ax, the having hence two salinity, The a time s c a l e can be ignored, motion that a value greater and I n a d d i t i o n , we than one; s p e e d s c a l e and must t h e r e f o r e assume a t t h i s p o i n t as t w o - d i m e n s i o n a l momentum i s p r i m a r i l y geostrophic. where and thickness. waves w i t h a R o s s b y number much l e s s U f / A x << 1 km. local coastal formulating ocean d e n s i t y , and i s described equations. motion described only nearshore u p p e r l a y e r of v a r y i n g d e n s i t y conservation than the l a y e r of c o n s t a n t upper l a y e r of information. considers l a y e r s ; a lower and in c o m p l e t e d e r i v a t i o n and that numerical a p h y s i c a l oceanography the model i s a v a i l a b l e i n [ 6 ] . i s a c o n d e n s a t i o n of The model is brief little the s e t of be of that f the > 10 want is, Coriolis much g r e a t e r Ax the km. than is suf f i c i e n t . Under these upper l a y e r assumptions, the are (1) -fv (2) fu = = (-1/p,)9p,/9x (-1/p,)9p,/9y ~ - Fx Fy momentum e q u a t i o n s for the 5 where: (u,v) = ( u ( x , y , t ) , v ( x , y , t ) ) = c o m p o n e n t s of p-, = pressure p, = density i n the of the upper l a y e r ; upper l a y e r ; ( F x , F y ) = c o m p o n e n t s of Since the lower assumption thickness and of (2) can the the be surface can i n t e g r a t e d over the (H Ap/p,) - 2 f U = -g/2 (H Ap/p,) at r e s t , height u p p e r l a y e r and (3) - f V = -g/2 (4) friction. l a y e r i s assumed t o be that velocity; x - 2 y is and much therefore be using less the than the neglected, (1) upper l a y e r t o yield Fx Fy where, (U,V) = where (U(x,y,t),V(x,y,t)) transport is = U = c o m p o n e n t s of velocity integrated transport, over the upper layer; H = H ( x , y , t ) = the p 2 Ap thickness = d e n s i t y of the = p - p, 2 of the upper l a y e r ; lower l a y e r ; = d e n s i t y d i f f e r e n c e between l a y e r s ; g = gravity; In the equations a b o v e , we have a d o p t e d the partial d e r i v a t i v e s ; they are thesis by 3 /3x s u b s c r i p t s of In (5) (6) the terms conservation (for + V -U U + example) transports, equations (HS ) i n d i c a t e d h e r e and as in notation for throughout this equation (1) or by independent v a r i a b l e . of H usual = + V^US are W e = W S e b the volume and salinity 6 where: S Sfc(x,Y/t), = M salinity S^ = a c o n s t a n t W V e K = the entrainment = the divergence Density equation a constant can of (7) where lower given Ap = p , a ( T ) a for the S(x,y,t) equations Here S(x,y,t) salinity salinity, using (3) through (6) (9) fU = - ( g a / 2 p , ) w i t h an + V .U= W lower to (10) H (11) (HS)^ + V -US = 0 - relate rewritten (H S) - kU X 2 y - S and can be 2 and is assumed encountered. f e (7) (H S) T, = S (x,y,t) between -fV = -(ga/2p,) I A (x,y,t) upper density , the layers (so to salinity, as kv € K have taken to the friction coefficient, produce transports. combining of temperature, temperatures and proportional friction components components of and i n t e g r a t e d Multiplying (Fx,Fy) velocity over (8) the by k (u,v), to with upper be k the layer and (9) to by f and gives (12) (f 2 (13) (f 2 + k ) U = -(ga/2 2 + k Substituting C = of r a n g e of K a function a p p r o x i m a t e l y by (8) we salinity; operator. > 0), t layer; velocity; function in upper S* Introducing difference layer be c o n s i d e r e d state is in the ) V = -(ga/2 2 (12) gak/(2p,(f +k )) 2 P l 2 and (13) ) P l ( f(H S) 2 ) + k(H S) ( -f(H S) 2 into ) 2 y x + k(H S) 2 ; < (10) y yields, ) with 7 (14) H Similarly, = W t + CV ( H S 2 t 2 substituting (HS\ ) . (12) and (13) i n t o (11) g i v e s = C (kSV (H S) + fJ(S,H S) + kVS-V(H S)) 2 2 2 2 or (15) S = C/H t ( k S V ( H S ) + f J ( S , H S ) + kVS V ( H S ) ) 2 2 2 where V i s the Laplacian operator, 2 i n terms of x and function field; H S currents nonlinear, governing Solving while y and represents 2 we e x p e c t causing (S/H)H V is the gradient distribution operator operator. Both differential f o r our model o f t h e (14) and equations coastal S(x,y,t), the salinity i n H S, 2 (15) a r e and a r e t h e circulation. (14) y i e l d s H ( x , y , t ) , t h e t h i c k n e s s of t h e upper s o l v i n g (15) g i v e s The o f t h e d e n s i t y mass to create gradients i n t h e model d o m a i n . equations t J i s the Jacobian the f r e s h water r u n - o f f parabolic partial 2 layer, defect, the d i f f e r e n c e b e t w e e n S^, a n d S . b The and parameters assigned values g = 9.8 m / s e c and c o n s t a n t s as f o l l o w s : 2 f = 1.1E-4 s e c " 1 k = 0.25E-4 s e c " a p, W =0.77 p p t " = 10 3 i n (14) a n d ( 1 5 ) a r e d e f i n e d 1 1 ppt = 6.0E-5 F v 2 e where F 2 = F r o u d e number = ( U v = surface v e l o c i t y 2 + V )/(gaSH )/p, = (U 2 2 + V 3 2 )''*/H 1 8 Further the d i s c u s s i o n of p a r a m e t e r s k, values transport, chosen (U,V), algebraically 2.2 The from Domain Of The domain ocean t o the can be found and hence the (12) and The Model of the and idealized for the for c o m p o n e n t s of can be found south, and b o u n d e d by on t h e e a s t by Strait coast Islands (see of J u a n de F u c a and Peninsula north I t i s through run-off enters coastal open the i s b r o k e n between t h e O l y m p i c of these waters. m o d e l , t h e d o m a i n has been t o a square, w i t h x d e s i g n a t i n g the e a s t - w e s t a x i s a x i s , as 2 . 3 Boundary by F i g u r e of t h e t o be p a r a l l e l 2. The square. and origin The to the y-axis. (12), requires land Conditions equation the g r a d i e n t illustrated southwest corner boundary i s t h e r e f o r e taken f o r the salinity on a l l b o u n d a r i e s be defect, zero; that that i s , 3S/3n = 0 where 9 /3n the v, i s an a r e a t h e p u r p o s e s of a n u m e r i c a l t h i s domain i s the (16) The velocity I s l a n d by Queen C h a r l o t t e Sound. north-south The [6]. t o t h e Queen C h a r l o t t e the passages that r i v e r For y coast V a n c o u v e r I s l a n d by Vancouver two The in justification (13). n o r t h , west and 1). and equations from the Olympic P e n i n s u l a Figure a, W no-flow i s the normal d e r i v a t i v e . condition boundary i s broken, across representing a This i s the boundary. fresh domain, S i s p r e s c r i b e d from h i s t o r i c a l water data. usual Where flow f o r m of the into land the 9 The It is boundary c o n d i t i o n s f o r H are a l i t t l e assumed more difficult. t h a t n e a r t h e open o c e a n b o u n d a r i e s , c h a n g e r a p i d l y and can therefore be approximated H does from those The actual v a l u e s w i t h i n t h e domain t h a t a r e near the boundary. method f o r f i n d i n g H on the numerical the these scheme has x-direction On + k(H S) x Expanding the d e r i v a t i v e terms + fH Sy boundary, be no flow (12), yields + 2kSHH 2 y the l a n d when = 0 2 y be p r e s e n t e d z e r o , as t h e r e c a n t h e r e f o r e , from e q u a t i o n 2 2fSHH will been e x p l a i n e d . t r a n s p o r t must be a c r o s s l a n d , and f(H S) boundaries not + kH S^ = 0 2 ;t or (17) Hy To find H = -(fHSy on the approximated (15). differential r e g i o n s of t h e The thickness constant kHS )/2fS x boundary, and a Sy first we note H that H can A be i s known f r o m t h e s o l u t i o n order initial of H i n y, w i t h t h e the v a l u e of given at value ordinary i n i t i a l value the river of for the discharge boundary. solution historical + x Conditions H ( x , y , 0 ) and initial is equation being 2 .4 I n i t i a l land numerically T h u s , (17) equation + 2kSH of (14) and S(x,y,0). river and values. field For requires i n i t i a l conditions S i n c e t h e model i s i n t e n d e d run-off salinity (15) data, data some will will be with average of e v e n t u a l l y be c o n s i d e r e d the p u r p o s e s of our values historic f o r use used i n i t i a l computer for both as runs, H(x,y,0) and 10 S(x,y,0). 2.5 How The M o d e l W i l l Equations specified, (14) and (15), with describe the effect of circulation. approximate data Be U s e d river The p u r p o s e of t h e model actual currents these to try values, t r a n s p o r t s ) at the inflow regions at Similarly, daily lighthouses series for thickness surface A along use as values data velocities (and were c a l c u l a t e d and • a surface salinity and hence weekly f o r t h e model readings, recorded t h e c o a s t , have been s m o o t h e d t o a w e e k l y salinity have been boundary values. approximated Upper layer from t r a n s p o r t s and salinities. period interpolated available of approximately weekly twenty t r a n s p o r t and s u r f a c e years of salinity actual of the study that l e n g t h of time, migration or records are t o be u s e d a s b o u n d a r y c o n d i t i o n s f o r t h e m o d e l . aim salmon closely discharge t i m e s e r i e s p r o d u c e d t o be u s e d a s b o u n d a r y v a l u e s equations. to o f W a s h i n g t o n S t a t e , B.C. surface as by u s i n g c o l l e c t e d To t h a t e n d , r i v e r h a s been c o l l e c t e d f o r t h e c o a s t s From is conditions r u n - o f f on c o a s t a l and s a l i n i t i e s i n boundary c o n d i t i o n s . Alaska. boundary The i s to approximate the c o a s t a l c i r c u l a t i o n f o r and t o compare t h e r e s u l t s w i t h data. From this i t is hoped the that known some i n f e r e n c e s c a n be drawn a b o u t t h e r e l a t i o n s h i p between t h e two. 11 Figure 1 - Domain of the Coastal Circulation Model Figure 2 - Domain of the Numerical Model 13 III. 3.1 CHAPTER 2 Notat ion To s o l v e the p a r t i a l (1) H t = CV (H S) + t = C/H 2 (2) S numerically, equations W 2 e ( k S V ( H S ) + f J ( S , H S ) + k V S V ( H S ) ) - (S/H)H_,. 2 the rectangular differential 2 domain grid. 2 is 2 discretized Initially, by overlaying t h i s domain w i l l be a a simplified v e r s i o n o f t h a t shown i n F i g u r e 2 by c o n s i d e r i n g i t t o have just a single discharge point. over the Ax The model e q u a t i o n s r e g i o n r e p r e s e n t e d by t h i s g r i d , w h i c h and are LAx Ay i n t h e x and specified = MAy = 600 dividing the specific time The grid by km. by domain i n N e q u a l (1) and the time or (2) w i l l discrete respectively. each finite finding numerical There (1) and and are of (1) difference the grid j=1,...,M, and discretized by s t e p s of s i z e A t . be approximated solutions For and A ease at any H(iAx,jAy,nAt) and of notation, the and (2) i n d i s c r e t e t e r m s , in some way. the Finite methods a r e u s e d most o f t e n f o r solutions to p a r t i a l a number of f i n i t e ( 2 ) , but on of nAt. p a r t i a l d e r i v a t i v e s must be a p p r o x i m a t e d element spacings is also s o l u t i o n s a r e s i m p l i f i e d t o H-- specify grid Coordinates time d i r e c t i o n i s r e p r e s e n t e d by S(iAx,jAy,nAt), To The time point valid ( i A x , j A y ) , w i t h i = 1 , . . . , L and s o l u t i o n s of approximate y directions. has are differential e l e m e n t methods f o r equations. approximating t h e y a r e f o r t h e most p a r t d i f f i c u l t i n general are computationally expensive to t o program solve. Finite 14 differences, costly to in addition solve, to being approximate w h i c h a r e more s u i t e d 3 .2 Finite Differences will use space d e r i v a t i v e s They will first and for and the first using second order x - d i r e c t i o n divided f o r the order first time second derivative. and d e r i v a t i v e s , and and Dj D f o r the x and ey Dj y derivatives. some f u n c t i o n D„X u i j less r e c t a n g u l a r domain. be d e f i n e d u s i n g l i n e a r o p e r a t o r s D^ the y - d i r e c t i o n For t o our centred differences order t o p r o g r a m and derivatives differences We easier = (u„ . u, u = u ( x , y , t ) , - u-„ . t| ( let )/2Ax and D ujj *x From = ( u Taylor ltt) K u = D X (i )/(Ax) expansion 2 of u((i+1)Ax,jAy,t) and that and 2 = D„ u + 0( (Ax) ) 2 X 2 2 + u;_ . i t i s e a s i l y seen 2 xx; 2UJJ u + 0((Ax) ) ox where 0 ( ( A x ) ) D - series u((i-1)Ax,jAy,t), u - linear operators approximations. the y - d i r e c t i o n Using the difference i s the o r d e r of the said Similarly, D derivatives to and oyf and in (1) in can and be (2). time D^ ahead some t=(n+l)At d e p e n d on v a l u e s a t l o w e r D^ and accurate are l i n e a r operators f o r , D| x , D^ solutions marching The second order constructed The in be term. are accurate to l i n e a r o p e r a t o r s D^ methods derivatives are remainder 2 and to Djy , finite approximate H and fashion; time 0((Ay) ). S are the levels. found solutions the by at 1 5 If H and previous S are found time l e v e l s , individually, t h e method i s said computationally very simple, e x p l i c i t using values from t o be e x p l i c i t . While methods a r e r e s t r i c t e d some maximum s i z e o f t i m e s t e p u s e d t o a d v a n c e i n t i m e . to Using a t i m e s t e p g r e a t e r t h a n t h i s maximum, w h i c h d e p e n d s on Ax a n d Ay, results in solutions instabilities severely Chapter by using restricts 1, that a grow e x p o n e n t i a l l y . time step less than t h e c h o i c e o f Ax a n d Ay. various assumptions used in A v o i d i n g such the maximum As was m e n t i o n e d i n deriving e q u a t i o n s p u t c o n s t r a i n t s on t h e s i z e o f Ax a n d Ay. method i s t h e r e f o r e n o t a p p r o p r i a t e f o r f i n d i n g the An model explicit solutions t o our equat i o n s . A finite dependent d i f f e r e n c e method i s i m p l i c i t variable interrelated. points Implicit system of l i n e a r taking at more at methods thus e q u a t i o n s a t each computing time the next A required time l e v e l . than implicit solution While explicit method e q u a t i o n s i s t h e C r a n k - N i c o l s o n method. defined level the of a naturally methods, on are implicit time step, f o r accuracy. well-known operators time require methods g e n e r a l l y do n o t p l a c e r e s t r i c t i o n s beyond t h a t i f the values of the above, for solving parabolic In terms of t h e the simple two-dimensional linear diffusion equat i o n u is t = KX U approximated + u yy by the Crank-Nicolson ( o r CN) method at ( i A x , j A y , ( n + l / 2 ) A t ) by - u|j = A t ( D £ ( + u ] )/2 + D/ { y ( u?j' + ujj ) / 2 ) 16 The difference second scheme i s c e n t r e d a b o u t t = n + l / 2 t o p r e s e r v e order accuracy equation, The solving finite (1) and (2) w i l l Before in the be limited two-level remainder t o t h e CN a p p l y i n g t h e CN difference t=n. of the their chapter which application to method. scheme t o e q u a t i o n s (1) and ( 2 ) , we t h e s p a t i a l d e r i v a t i v e s a r e t o be e x p a n d e d . equations have variables which derivatives must be of products expressed as an e x a m p l e t h e s e c o n d ( H S ) ^ , in equation (1). 2 u s u a l way a d i f f e r e n c e e q u a t i o n s and must d e c i d e how Consider is f o r t=n+1 f r o m v a l u e s a t discussion concerns of the as d i f f e r e n c e The dependent equations. order x - d i r e c t i o n d e r i v a t i v e Applying the product rule in , the yields (3) 2 H S H while i n t i m e , and the + 2H S + 4HH S 2 K)C X considering X (H S) A + H S^ 2 as a p r o d u c t 2 of f u n c t i o n s (HS) and H gives (4) (HS)H X;< + 2(HS) H^ + H(HS) X Upon d i s c r e t i z a t i o n , except for fortuitous s o l u t i o n s of both d e r i v a t i v e s are We that expanded runs (1) and (2) w i l l functions H and clearly S. unequal, Obviously, depend on w h i c h (4) functions. above equation (1) f o r H, multiplying was (as e x p l a i n e d i n Chapter Expanding of are way the the (H S) 2 expanded. separate as expressions values are s o l v i n g equation those single, these X X That H the and i t seems should reasonable be c o n s i d e r e d derivatives should c o n f i r m e d by e x p e r i m e n t a l 4). (1) t h e r e f o r e y i e l d s as be computer 17 (5) H The = C ( ( H S ) V H +2(HS) H ;< H V ( H S ) t e r m may be 2 is, V (HS) is 2 ( [ 1 2 ] , p. + 2(HS) H 2 t 195) / t r e a t e d as considered the the t o assume t h a t dependent v a r i a b l e are ;t known. (5) s i n c e we than are that It is usual the K + (HS) H )) + y y f, e for S also contains s o l v i n g f o r S, outlined f o u n d t o work that +W . 2 equation terms i n ( 5 ) ; becomes X The e inhomogeneous terms i n v o l v i n g 2 = HV (HS) 2 c o e f f i c i e n t of H. (6) R t = C ( ( H S ) V H + 2 ( ( H S ) H where f , other +HV (HS)) + W . y their f o r H. best, and d e r i v a t i v e s of expansion Treating thus, for H must (H S), 2 be different as a s i n g l e e n t i t y 2 (H ^)^, the 2 but was appropriate expansion i s H S 2 + S(H ) 2 X X X ; < Expanding equation (7) = C/H + 2S (H ) 2 X (2) ( X gives ( k S ( H V S + 2VH -VS + S V H ) + 2 2 2 2 2 fSJ(S,H ) + k (SVS-V(H )+H (VS) ) 2 Just as as f o r H, t e r m s not t = C/H 2 (7) ( k S ( H V S + 2VH * VS) 2 2 2 = (CkS/H) ( S V H ) 2 2 W h i l e t h e method of determined, obviously are the 2 nonlinear. nonlinear iterative of 2 2 2 <; expanding the derivatives the d e r i v a t i v e s t h a t the dependent v a r i a b l e s ; Discretization of such both equations systems of e q u a t i o n s w h i c h would r e q u i r e scheme f -(S/H)H . c o e f f i c i e n t s of functions treated + fSJ(S,H ) + 2 2 t becomes k(SVS«V(H )+H (VS) ) + where f (S/H)H 2 i n v o l v i n g d e r i v a t i v e s of S a r e f o r c i n g f u n c t i o n s , and (8) S 2 has been result are equations produces some k i n d f o r s o l u t i o n , s u c h as Newton's method. of Since 18 iterative methods for computationally linearizing very the having 2, equations expensive, we are in approximate general further by equations. A l t h o u g h we Chapter nonlinear leave the variable the result computational of linearizing coefficients. details (6) and The e q u a t i o n s until (8) can later are in equations be simplified to and s -t = where a , of linear i xx s b^ , K the b + b S b2Syy + i = 1,..,4, operators 3 are + b^Sy + x functions defined of earlier, f 2 x, y and t . the CN In terms discretization i s (9) Hf* - H-5 a D = At(a^H^+ 2 2 H?; - a D ^ H ^ Y y + + a D a 3 H^ H f,) y o y and (10) = A t ( b , D i S ^ + S?f ~ S% b 2 D J y S i ^ + b«D It each should derivative values both of equations ([12], 26) of the here the that equations and u s i n g are the the evaluating constants diffusive. difference scheme central terms if diffusion insignificant, terms. the o y x S^+ ( f ). 2 coefficients for defined in Chapter 1, in is important because it is are coefficients While the p r i m a r y moving poor are large is well known approximations convective force of representative This differences their 3 the shows t h a t , strongly that convective those not p. of H and S, determining to be n o t e d b EL S ^ + relative terms p r o d u c e d by to are the 19 diffusion terms. 3.3 Consistency, S t a b i l i t y Finite difference approximations show that the approximation we to Convergence equations differential finite to And naturally equations. difference the exact are I t i s necessary solution solution. for A x , A t -> 0; that the d i f f e r e n c e scheme w i l l give a same s o l u t i o n ) ; and for the solution w i l l exact approach initial and standard theory does the (assuming s o l u t i o n as A x , A t -> and and finite the difference 0. These convergence. convergence of coefficients, not c o n s i d e r the e f f e c t of e i t h e r boundary c o n d i t i o n s model problems the chapter and of t h a t the f i n i t e stability 'good' differential v a l u e p r o b l e m s g e n e r a l l y assumes c o n s t a n t or v a r i a b l e c o e f f i c i e n t s . the for to t h i s means t h a t solution p r o p e r t i e s are c a l l e d consistency, s t a b i l i t y The a the solution bounded the exact is Briefly, r e q u i r e the d i f f e r e n c e equation to tend to equation only We will assume f o r meet t h e s e c r i t e r i a , the and moment that consider later t h e e f f e c t on o u r a n a l y s i s o f v a r i a b l e in coefficients boundary c o n d i t i o n s . To simplify scheme, we t h e p r e s e n t a t i o n of t h e s e c o n c e p t s c o n s i d e r t h e one equations. dimensional version What f o l l o w s i s e a s i l y e x t e n d e d t o two For u = u ( x , t ) , l e t (11) w i t h a and u t = au ^ x + bu„. + b c o n s t a n t , and u(x,0) = g ( x ) . of f(x,t) with i n i t i a l condition f o r the the CN model dimensions. 20 Much of and the analysis that definitions the CN scheme t o of K r i e s s and (11) (1 - a A t D ^ follows Oliger initial where v ( x , t ) any - bAtD )v(x,(n+1)At) = i s the the matrix C ! X s o l u t i o n t o the discretization notation; )v(x nAt) + and (11) can finite of v(x,At). be The particular 2 c o n s t a n t c and difference expressed for a l l s u f f i c i e n t l y can is accurate u(x,t), if ||» || We i s the have seen a p p r o x i m a t e 9/9x CN 1 small At, a p p r o x i m a t i o n of h(x)=g(x), that the series error i n our truncation in the 9 /9x 2 9/9t 2 with is a interval Ax + Atl* ) For operators, accuracy 0 ( ( A x ) ) , 2 i s accurate to 0 ( ( A t ) ) . initial the linear condition and is scheme, c a u s e d by the D,* and Dj^ , and that the 2 Assuming that known exactly, the (2,2). a p p r o x i m a t i o n s t o the error. there + At?*" ) l scheme i s a c c u r a t e of o r d e r The order norm. 2 that and of every f i n i t e 0 - g (x) | | < c ( A x t s i m p l i f y these as | |Q,u(x,t+At) - Q u ( x , t ) - A t f | | < C ( t ) ( A x V | |h(x) linear Atf(x,nAt) scheme solution equation. produces We a f u n c t i o n C ( t ) , b o u n d e d on [0,NAt], such t h a t difference (11) 0 ( q , q ) f o r the where + Atf(x,nAt) f Q,v(x,(n+1)At) = Q v(x,nAt) + Definition ( Applying h(x), x, (12) 29). (1 + a A t D ^ w c o m b i n a t i o n s of v ( x , ( n + l ) A t ) using notation condition v(x,0) = For ( [ 5 ] , p. the yields bAtD with i s b a s e d on CN truncating derivatives, is scheme, t h e truncation the Taylor called error the i s At 21 times local the accuracy, or error, the error 0(At +AtAx ). 3 2 i n the s o l u t i o n Note t h a t t h i s from one t i m e i s the step to the next. If the d i f f e r e n c e equation converges equation -being consistent. from the approximated, F o r t h e CN definition bounded, and q i , q the scheme, of accuracy. the is differential said to follows be directly S i n c e C ( t ) i s assumed t o be > 0, 2 a s A t , A x -> 0 1 t h e scheme i s c o n s i s t e n t . While the difference solution s o l u t i o n w i t h t h e order of accuracy depends (Ax) 2 on terms 2 these several steps solution. in of c o u r s e are First, taken the equal to nondimensionalize when the At(Ax) . 2 Stability error expansion. To m i n i m i z e finding In the e r r o r , a numerical step; thus (At) will be i t i s conventional to 3 This requires the equations t o The t i m e scale now depends on used. requires equation, v ( x , t ) , exact t h e t i m e s t e p i s c h o s e n t o be Second, the equations. resolution series actually space s c a l e d so t h a t 0 < x , y < 1 . the s p a t i a l Taylor of the d e r i v a t i v e s m u l t i p l y i n g the a r e unknown. size a b o u t t h e same s i z e a s roughly the approximates g i v e n above, t h e exact the s i z e of the p a r t i a l and ( A t ) practice, be method consistency C ( t ) ( (Ax)**' + ( A t ) * " ) -> 0 and to that the solution to the d i f f e r e n c e a t t = n A t , i s bounded a s A t , A x -> 0. Def i n i t i o n The d i f f e r e n c e a p p r o x i m a t i o n ( 1 1 ) is stable s e q u e n c e A t , Ax -> 0 i f t h e r e a r e c o n s t a n t s a , Kj s u c h s t ^ 0, t h e e s t i m a t e for a that f o r 22 ||v(x,t)|| = | |S(t,0)v(x,0) | | < K exp(a t) | | v(x,0) | | s s holds. Here S ( t , 0 ) i s a s o l u t i o n operator with S(0,0)=I, S ( t , 0 ) = S ( t , t ) S ( t ,0); t >t,>0. 2 For t=nAt, (13) and (12) 2 1 l 2 is Q,v(x,(n+1)At) = Q v(x,nAt) + Af(x,nAt) 0 S ( n A t , 0 ) can easily be found. F o r Q = Q7'Q , 0 v(x,(n+1)At) = Qv(x,nAt) + Qf Atf(x,nAt) 1 = Q(Qv(x,(n-1)At)+Q f A t f ( x , ( n - 1 ) A t ) 1 + Q ~ A t f (x , n A t ) 1 It is condition evident in contribution should the from the definition stand For factor that is. i s not from greater one time w, we t h a n one. step to to the That i s , the bound we However, let the h a v e t o show t h a t hold. of must show t h a t That the is, v(x,t) exists. the a m p l i f i c a t i o n stability requires n e x t t h e a m p l i t u d e of solutions. is von include for simplicity, increase, except to the stability A t ^ S ( n A t , i A t ) f (x , i A t ) . f r e q u e n c y does not This the should that a Fourier transform frequency, of I t i s c l e a r t h a t we i s bounded f o r s t a b i l i t y every side above t h e b o u n d , and as Assume f i r s t right f o r c i n g f u n c t i o n , f. i n c l u d e u p p e r l i m i t s on increases -1 the definition t h i s only | |Q| I** that allow truly every increasing Neumann n e c e s s a r y c o n d i t i o n for 23 stability ([12], Applying transform ([9], p. D 70). Fourier of the 37), they = i PX p. transforms linear to (13) operators D requires and 0)l D . only the As shown 2 e in are sin(/3/Ax) (14) D. , = -4 2 and G , the sin (/3/2)/(Ax) , 2 2 amplification 1 - factor, is (from [8], 2asin 0 + i (b/2) ( a A t / a ) sin(2/3) 1 + 2asin 0 + i(b/2)(aAt/a) sin(20) 2 p. 196), + Atcf G = 2 with a = a A t / ( A x ) equation (1 in the and 2 j3 = k A x . + 2asin /3) 2 2 sin (20) it greater 2 denominator; - + (b /4)aAt/a is modulus, (15) the to an 1 s i n (20) 1 f o r any /3. So + Atcf) i b / 2 ( a A t / a ) sin2/3 ) factor can be g i v e n by 3 2 2 1 + O(At), h (/3), 2 h (|3) and h ( j 3 ) 3 a exponential solution b e c a u s e of satisfies the the it is difference schemes for shows that ([12], scalar as the are bounded f u n c t i o n s . that may e x i s t function. condition, known sufficient growth forcing von Neumann stability; This 2 than (l+h (/3)At+h d3)(At) +h3(/3)(At) +h«(0)(At)M term a l l o w s is + 2asin /3 - amplification |G| < h,(|3), condition clearly 2 = O(At) leads 2 2asin /3 + i ( b / 2 ) ( a A t / a ) (1 where G with G < (1 In Simplifying p. in the true The i n e q u a l i t y above which 72) The that equations, is necessary for a l l the for two-level von Neumann well. solution from one time level to the 24 next is bounded. ||S(t,0)v(x,0)|| bounded, and stability definition. values of the of At unconditionally If and the CN scheme d i f f e r e n c e equations tends Assume a c c u r a c y and that stability ||u(x,t)-v(x,t)|| Since i t was that stability At, -> Ax 0, and was equation |G| i s bounded scheme is the for said s o l u t i o n s to the zero 0. as At, convergent. fixed At and are valid. any to be Ax < e x p ( a t ) (Ax*' +At* scheme, T h e r e f o r e the convergence Ax -> 0, a From [ 5 ] , p. 31, we the 1 5 CN differential d e f i n i t i o n s for Then f o r a l l t = n A t , i s unconditional, there n > 1, )c . that (q,,q )=(2,2) and i s no restriction the 2 r i g h t side tends to zero on as At, i s proven. mentioned, the with constant equations c l e a r l y do in stability One terms evaluated their a n a l y s i s above assumes a coefficients, not meet. variable coefficients, errors in of requirements for shown, f o r t h e -> be Variable Coefficients As by the to scheme i s s a i d t o be Theorem 3.4 and also stable. following: Ax Ax, therefore the I t i s c l e a r that and have t h e way meets d i f f e r e n c e between the difference must 'local condition accuracy a s s u m i n g of c o u r s e t h a t it is the is model unaffected there are no quite possible that the change. method o f e n s u r i n g of W h i l e the calculation, c o n d i t i o n may a differential stability stability'. at each g r i d p o i n t i n the i s to look That i s , the domain and at the problem c o e f f i c i e n t s are tested separately 25 for stability. usual that locally W h i l e t h i s does not instablities ( [ 1 2 ] , p. in 91). guarantee convergence, i t i s the approximate However, t o a p p l y t o a l a r g e number o f p o i n t s p r e s e n t s practical way to test for in 3.5 c a r r i e d out Chapter stability Of Boundary stability, we must p r o d u c e d by the CN check that the the m a t r i x Q any to a by are is to the This r e s u l t s are the boundary system of initial given conditions linear value z with on |z| > assume that s u c h an In p a r t i c u l a r , we 1, t h e a l l boundaries. at the of eigenvalues so o n l y consider the the of the eigenvalue to z is this leads and the (12), this for boundary c o n d i t i o n s no-flow boundary be with H for S condition of H and corresponding equation add conditions no-flow c o n d i t i o n w i l l sample not equations. Assuming smoothly changing f u n c t i o n s change the do and want t o show t h a t boundary inflow regions, equations problem, show t h a t t h e model The on 1. boundary l a y e r p r o b l e m s ) , p r e s c r i b e d p r o b l e m , and We the boundary c o n d i t i o n s prescribed not of b o u n d a r y c o n d i t i o n , and s p e c i f i e d i n Chapter elsewhere. do examine eigenvalues the prescribed t h u s no problems. i n t r o d u c t i o n of b o u n d a r y c o n d i t i o n s o c c u r s f o r the are effect contradiction. As sample Conditions scheme f o r t h e strategy introduced One of v a r i a b l e c o e f f i c i e n t f o r s e v e r a l p r o b l e m s , and determine the The Neumann's method 3. Stability To von start obvious d i f f i c u l t i e s . p r o b l e m s i s t o p e r f o r m c o m p u t e r r u n s of was solution S (and conditions initial value analysed. no forcing 26 function, and problems, 0 ^ x < 1. stability stability of f o r both It the left boundaries the value a t u ( 0 , t ) w i l l is boundary that for parabolic condition implies [ 1 8 ] , and so o n l y t h e e f f e c t o f be c o n s i d e r e d . second order a c c u r a t e approximation the l e f t known In terms of u ( x , t ) , the t o the no-flow condition boundary i s u(0,nAt) = (4/3)u(Ax,nAt) - (1/3)u(2Ax,nAt) We have s e e n t h a t t h e a m p l i f i c a t i o n m a t r i x , Q = Q p Q o r have at eigenvalues must n o t g r e a t e r t h a n one f o r v ( x , ( n + l ) A t ) t o be a bounded s o l u t i o n . L e t g = (g*,...,g«) be an e i g e n v e c t o r of Q, so t h a t Qg = z g or z Q i 9 = Qo9 with z an eigenvalue. their coefficients, TL c z ( 16) £ From [ 1 8 ] , and Q 0 c a n be e x p r e s s e d = X.a g; ; - a )gj_ = 0 . t (16) or is a linear difference c o n s t a n t c o e f f i c i e n t s , and t h e s o l u t i o n 9 where (16) = the where equation £aj(T.(z)) Tj ( z ) a r e t h e r o o t s o f ( 1 6 ) . . S u b s t i t u t i n g t h i s 2 - 2/3T B= with t a k e s t h e form, gives (17) T a s sums o f s o i n t e r m s o f g a n d z, g; jL(zCj Q, + T/TJ = 0 1/(1+Ax)+((z-1)/(z+1))/(X(1+Ax)) r = X/2 - XAx/2 T? = X/2 + XAx/2 into 27 and X = For | z | > 1, with At/(Ax) |T,(z)| < (18) and, qi to (17) 1, an e Substituting (19) The = be also T h u s , we (4/3)g, (19). c o n t r a d i c t i n g our t h a n one. and T 2 Thus z, w i t h 2 the no-flow that . (18) gives + (l'/3)T^z) ) = 0 p r o d u c e an The finding must a l s o s a t i s f y require into 1 - (4/3)T,(z) r o o t of gj - d/3)g |z|>1, an e i g e n v a l u e , a T, So eigenvector, boundary c o n d i t i o n w i l l z, w i t h roots (z) this expression a,( shown t o have 2 boundary c o n d i t i o n . g can | T ( z ) | > 1. = a,T, be 2 unstable solution i f , for the corresponding r o o t s of t h a t one r o o t of (19), however, are of the |z|>1, i s n o t boundary, approximated to second order an r o o t s of (17) eigenvalue; accuracy, (17) is 1 and 3, is the does not less no-flow affect stability. 3.6 Nonlinear To it avoid having i s necessary terms. Linear but the the true To the the t o s o l v e our to l i n e a r i z e equations linearization nonlinear equations the c o e f f i c i e n t s of the a r e , of c o u r s e , directly, derivative much e a s i e r t o p r o c e d u r e y i e l d s o n l y an solve, approximation to of the we coefficients. examine note f i r s t of Terms this e f f e c t on the e q u a t i o n s the accuracy are q u a s i l i n e a r , s i n c e the second d e r i v a t i v e terms i n v o l v e only function i t s e l f . Although t h e CN first system, coefficients d e r i v a t i v e s or d i f f e r e n c e scheme can be 28 modified to truncation produce error predictor-corrector estimate is formulae 0(AtAx +A\? -) 2 [ 8 ] , the rather fx than 0(AtAx +At ) w h i c h t h e CN scheme ( f o r l i n e a r e q u a t i o n s ) y i e l d s . Lees has 2 3 [8] coefficients shown that by using extrapolated estimating values from the the two time steps, the t r u c a t i o n e r r o r remains 0 ( A t A x + A t ) . CN i s centred 2 scheme about t=n+l/2, t e r m must be c a l c u l a t e d a t t h a t terms in the time. coefficients quasilinear previous Since the 3 t h e c o e f f i c i e n t s of each In p a r t i c u l a r , i f the of the second order HS partial d e r i v a t i v e s o f ( 9 ) a r e a p p r o x i m a t e d by US"* ' - = 1.5 H S 1 7 an e r r o r o f 0 ( ( A t ) ) i s added t o t o t h e t r u n c a t i o n 3 thus remains We involve error, which 0(AtAx +At ). 2 see that 3(HS)/3x same way. - 0.5 HS*",' w 3 the coefficients and 3(HS)/3y F o r (HS) a t p o i n t o f 3H/3x a n d 3H/3y, w h i c h t e r m s , c a n be a p p r o x i m a t e d (i,j), i n the the d e r i v a t i v e a t time n+1/2 c a n be a p p r o x i m a t e d by (HS) X = ( (HS)j^ - (HS)^- ) / (2Ax) (HS) y = ( (HS) ;-„ - (HS) ) / (2Ay) c l i > ( and where a l l t h e (HS) v a l u e s a r e c a l c u l a t e d a t t = n + l / 2 a s s p e c i f i e d above. This method c a n a l s o be a p p l i e d One of c a n be t r e a t e d a s t h e c o e f f i c i e n t o f t h e o t h e r , and (10) 2 x the same time of n+1/2 each the ( [ 2 ] , p. 1 4 1 ) . has q u a s i l i n e a r c o e f f i c i e n t s which a r e f u n c t i o n s of S H , H , H , Hy, S 2 to S. terms derivatives extrapolating for 2 in Equation equation y t h a t appear a p p r o x i m a t e d by the t o t h e S% a n d S fashion x and Sy. A l l of these as o u t l i n e d above, are estimated in without a l t e r i n g the order 29 of the t r u n c a t i o n e r r o r . The e f f e c t on procedure [19], stability of applying this to a l l the terms i n each e q u a t i o n however, has shown that for extrapolation i s not an known. equation involving n o n l i n e a r terms i n the f i r s t s p a t i a l d e r i v a t i v e of the and the f u n c t i o n i t s e l f , the s t a b i l i t y extrapolation of depends on t h e r e l a t i v e t h e d i f f u s i v e and c o n v e c t i v e t e r m s . not highly convective, conditionally diffusive, While we that, stable. i t for seems l i k e l y no computer chosen To coupled, solve space the f o r c i n g to steps, t h e model e q u a t i o n s , we independent While be w i l l be to are be highly maintained. the results i n Chapter this are 4) of show linearization stable. note from this s l o w l y from one f i r s t that they variables are are present S simultaneously, in scheme the only r a t h e r than a p p r o x i m a t e l y S i n c e t h e y a r e n o n l i n e a r as w e l l , one i n any case. are of Little approach. I n s t e a d , assume, as was varies found the C r a n k - N i c o l s o n d i f f e r e n c e solved d i r e c t l y functions. gained coefficients equations e i t h e r H o r S w o u l d have t o be a p p r o x i m a t e d w o u l d be with Equations equations. terms scheme was a p p e a r t o be c o n d i t i o n a l l y c o u l d be a p p l i e d t o s o l v e f o r H and new scheme analysis, time function, For equations which (described f u l l y and CN s i z e of the model stability i n t h a t t h e two two CN the that s t a b i l i t y runs produces r e s u l t s which 3 . 7 S o l v i n g The the Since both can p r o v i d e preliminary the of Varah time mentioned level in Chapter t o the next. 1, that If equation S (9) 30 is solved S i s i n f a c t s l o w l y v a r y i n g , t h e e r r o r of approximation small first, HS h a s t o be e s t i m a t e d (at l e a s t smaller t h a n t h a t c a u s e d by solving equation addition t o the error i s0 ( ( A t ) ) , small Once H ,wVa (10) would has w i l l be estimating H , as 2 As we have s e e n , t h e which, will only i n the o v e r a l l e r r o r of approximating H If cause a H. been c a l c u l a t e d , t h e a v e r a g e o f H^'and H* , -, c a n u s e d i n e v a l u a t i n g t h e c o e f f i c i e n t s o f e q u a t i o n ( 1 0 ) , rather than an extrapolated p o s s i b l e t o repeat then S**' Each to step the process, calculate may be predictor-corrector interaction one require). 3 increase (by e x t r a p o l a t i o n ) . using t h e new v a l u e s values A+ o f H* M method. t o the next. n e x t c h a n g e by o n l y for S ' , i ti s We an can In t iteration lessen practice, and +v new a p p r o x i m a t i o n s f o r H "'' a n d S considered iterating until changing Solving n+Vl of -. a t h e e f f e c t of t h e o f H a n d S by a l t e r n a t i n g t h e o r d e r iteration continue value. o f s o l u t i o n from i t i s possible t h e s o l u t i o n s f r o m one i t e r a t i o n some a r b i t r a r y ( s m a l l ) amount. o f H a n d S, t h e s o l u t i o n s w i l l to to the F o r smoothly converge, i ti s hoped, i n o n l y a c o u p l e o f i t e r a t i o n s . Applying produces t h e CN d i f f e r e n c e scheme t o two s y s t e m s o f l i n e a r each time s t e p . equations t h e systems t o be s o l v e d a g a i n . the m a j o r i t y of computation time w i l l equations, better particularly accuracy. The a s we i n c r e a s e consequences equations t h a t must be s o l v e d a t I f some method o f i t e r a t i o n change t h e c o e f f i c i e n t s o f t h e m a t r i c e s , require t h e model i s used, which each iteration will will I t seems l i k e l y that be s p e n t i n s o l v i n g linear the resolution to obtain are discussed i n t h e next 3 1 chapter. 32 IV. We now solve, c o n s i d e r the the need. The of C h a p t e r and S??' effort, methods of of 2 d e p e n d on (i-1,j), the l i n e a r i n the vector their solution and of (i,j-l) as w e l l as banded both matrices are with half elimination former The l a t t e r can be c o n s i d e r a b l y f a s t e r . converge, and r e q u i r e , d e p e n d s on are the w i t h x, result y and well-conditioned that they The H^ linear row The with conventional the For resulting where L i s t h e our equations, be any solved using either number of the advantage of r e l i a b i l i t y , thus the The r a t e at which number of iterative Noting but the there will is matrices can be. (We being the will remain n o n s i n g u l a r . ) T h i s e l i m i n a t i o n w o u l d a p p e a r t o be we the best method. they For an o f t e n converge our to n o n l i n e a r terms, little the iterative iterations t h a t t h e c o e f f i c i e n t s of of a p p r o x i m a t i o n s t, methods a iterative t h e c o n d i t i o n number o f t h e m a t r i x . matrix i f at a l l . has or methods. slowly, the (10) nonsymmetric. v a r i a n t of G a u s s i a n ill-conditioned would (9) and t h e way b a n d w i d t h L+1, banded, nonsymmetric m a t r i x can methods the s p a c i n g away i n e a c h (i,j) . number of g r i d p o i n t s i n t h e x - d i r e c t i o n . A each unknowns. grid must that i s , to solve for ( i ,j ) i n v o l v e s ( i + 1,j), ( i , j+ l ) , being we o p e r a t o r s and a p p r o a c h i s t o o r d e r t h e g r i d p o i n t s by matrix matrices for equations i n v o l v e o n l y g r i d p o i n t s one spatial direction; the i n t e r m s of o p e r a t i o n c o u n t s , s t r u c t u r e of t h e m a t r i c e s are ordered operators structure possible computational CHAPTER 3 say matrix and vary about how assume o f case, course- Gaussian 33 To calculate t h e amount of c o m p u t a t i o n a l we n e e d o n l y l o o k a t t h e o r d e r a n d matrices. If (so Ax=Ay), Although L i s t h e number the order is L the effort required, half-bandwidth of g r i d p o i n t s i n each and 2 the matrices i n i t i a l l y the i n " these these elements matrix, zeros. must be having order considered m u l t i p l i c a t i o n s ) necessary 0(L"). form be of solved alternating the next as nonzeros. i tto schemes e x i s t only direction o f work, For such a ( t h e number o f triangular form is f o r equations with the t h e model e q u a t i o n s which in elimination t h e amount o p e r a t i o n count t o reduce However, d i f f e r e n c e same can the 2 L+1. have z e r o s b e t w e e n t h e e l e m e n t s When c a l c u l a t i n g L , the direction half-bandwidth on t h e edge o f t h e band a n d t h e d i a g o n a l s , G a u s s i a n "fills of produce m a t r i c e s 0 ( L ) operations. They 2 are which called i m p l i c i t methods, and a r e t h e s u b j e c t of section. 4. 1 I n t r o d u c t i o n To ADI Alternating direction implicit, or ADI, methods p e r t u r b a t i o n s of d i f f e r e n c e methods f o r m u l a t e d d i r e c t l y difference often approximations arise by adding multi-dimensional affect the accuracy CN scheme (1 - XD the d e r i v a t i o n for u 2 0 x scheme ADI from t h e m e t h o d s most operators in a way to a t h a t does not o f t h e scheme b u t a l l o w s an e a s i e r s o l u t i o n . t = u^ operators i n t r o d u c e d i n Chapter (1) derivatives. difference difference To i l l u s t r a t e the to are - XD^)u*" of a + u y y particular , i s , using 2, = (1 + X D 2 X + XD 2 y )u* ADI the method, linear 34 where X=At/2. If the operator d i f f e r e n c e o p e r a t o r on b o t h (2) (1 - XD^ - XDiy X D^ D 2 sides, + XD 2 (1) becomes 2 C (1 Factoring (3) If we ( 2 ) i n t e r m s o f D£ , D| x (1 - XD 2 X introduce )(1 - X D 2 - XD (4) 2 C X 2 y 2 X D )u* . 2 y level b e t w e e n n a n d n+1, namely n ) u*** . u*** i s a dummy v a r i a b l e The method of i n the sense t h a t t h e v a l u e t o t h e s o l u t i o n a t any p a r t i c u l a r splitting for splitting this )u 2 ) u**» = (1 + X D Peaceman-Rachford formula (4), 2 2 i n t o two e q u a t i o n s , 1 methods + X D« )u**> = (1 + XD,, ) ( 1 + XD , ) u * . U" ** i s n o t an a p p r o x i m a t i o n time. 2 e y V variable of + XD gives y )u*** = (1 + XD 2 ( 1 - XD The X an i n t e r m e d i a t e t i m e say n+*, (3) c a n be ' s p l i t ' (1 + XD| 2 y to the D„ )u^' = 2 < i s added 2 e y used in i s called ( [ 9 ] , p. 6 0 ) . A l t h o u g h a two-dimensional formula w i l l (4) hereafter be the t h e r e a r e many parabolic equation referred to as like the ADI scheme. The level either advantage of (4) i s t h a t t h e s o l u t i o n a t t h e next ( e i t h e r u**' o r u**') , i s n the x- or y-axis. involves p o i n t s immediately that i s , the matrix found along t o the point i s tridiagonal. system c o m p a r e d t o t h a t f o r t h e CN scheme. 2 being of equations 2 by L 2 is only effort small system, t h e multiplications, 0 ( 1 / ) r e q u i r e d by t h e CN scheme. This c o n s i d e r a b l e savings i n computation to found; The c o m p u t a t i o n a l F o r an L t r i d i a g o n a l matrix requires only 0(L ) to parallel The m a t r i x t o be s o l v e d t h u s adjacent required to solve a tridiagonal lines time time w i l l compared n o t come 35 at the expense method. of As we w i l l accuracy analyzed error, stable, i n C h a p t e r 2. error of for a of the d i f f e r e n c e 0((At) +At(Ax) ) 3 and 2 is sample problem s i m i l a r t o t h a t Note t h a t we h a v e assumed t h a t this will stability show, t h e P e a c e m a n - R a c h f o r d scheme, l i k e t h e CN scheme, h a s a t r u n c a t i o n unconditionally or i n c a l c u l a t i n g the Ax=Ay. truncation For the r e s t of the t h e s i s , be t h e c a s e . 4.2 A c c u r a c y And Stability To c a l c u l a t e t h e a c c u r a c y consider an equation of the ADI scheme, we o f t h e same f o r m a s o u r model first equations, namely, f o r u = u ( x , y , t ) , (5) u where = au t a,b,c using X x + and d + cu bUyy + du x , v are constant c o e f f i c i e n t s . t h e ADI scheme i n t h e same manner a s a b o v e , (6) (1 - X a D - XcD„ 2 x x )(1 - XbDj (1 + X a D where X i s a s a b o v e . (1 - X b D 2 2 x Splitting - XdD* v )u M + To 1 + XcD e #Jf ) ( 1 + XbD 2 - XcD e;t calculate the error, Taylor produce D^ and 2 o y y + XdD j ) u* , 0 * = (1 + X a D + X c D ^ ) u* 2 )u*«" = (1 + XbD 2 e y we s e e f i r s t (X ab)D^Dj^ . 2 * (\ cd)D 2 o x D + XdD cy that the operators , ey )u* * + (X bc)D D 2 o x c a n be a d d e d d i r e c t l y the to the t r u n c a t i o n S i n c e each o p e r a t o r errors can easily 2 < y , T h e s e t e r m s do n o t r e s u l t f r o m s e r i e s e x p a n s i o n of t h e d e r i v a t i v e s , and t h e e r r o r s CN scheme. u"**' , 2 ( 6 ) we g e t a d d e d t o t h e CN scheme f o r ( 5 ) a r e (X ad)D gives * (1 - X a D ^ (5) - XdD y ) u ^ = y (7) the Factoring i s applied error they found f o r to both u and H be f o u n d by e x p a n d i n g u ** 1 1 in a 36 Taylor we series. Examining j u s t one of t h e s e t e r m s , (X cd)D D 2 , have (X cd)D D 2 Cancelling u 2 on b o t h s i d e s , truncation the (u + A t u + 0 ( ( A t ) ) = ( X c d ) D error error we see t h a t of t h i s term so 3 as t h e CN the t r u n c a t i o n is error a s shown f o r CN Stability w i t h t h e CN CN a p p l i e d u addition The 3 to the contribution similarly f o r ADI scheme d i s c u s s e d e a r l i e r ; o v e r a l l a c c u r a c y unchanged. the is 0((At) ). from the o t h e r o p e r a t o r s 0((At) ), D 2 found to to be i s o f t h e same o r d e r the added terms Consistency follows leave the i n t h e same way i n the l a s t c h a p t e r . f o r t h e ADI scheme. scheme c a n be s e e n by a g a i n Looking at the a m p l i f i c a t i o n t o e q u a t i o n ( 5 ) , we comparing f a c t o r , G, for have 1-2(a,S,+a S )+i((b/2)(a,At/a)S +(d/2)(a At/d)S„) 2 2 3 2 (8) G = 1+2(a S +a S )+i((b/2)(a,At/a)S +(d/2)U At/d)S ) 1 where S , = s i n 2 / 3 , a 1 2 S =sin 7j, 2 2 = b A t / ( A y ) , 0 = k,Ax 2 2 It is easily 2 3 S = sin2/3, the a m p l i f i c a t i o n additional terms real valued sin(2/3) terms and denominator and unchanged. Thus, bounded by one (5). Using that t h e ADI the But and ADI for f o r ADI, we products 196) these . At,Ax->0. note In that the transformed to sin2(0), of sin (7?), 2 terms a r e added t o b o t h the of (7); from Chapter the factor i s unconditionally g u a r a n t e e d by t h e a c c u r a c y and p. < 1, amplification the r e s u l t s 2 scheme a r e F o u r i e r numerator the |G| factor involving sin(2r?). ([12], 2 a a, = a A t / ( A x ) , a and n = k A y calculating of S = s i n 2 r y and 3 shown 2 inequality for stable ADI for 2, c o n v e r g e n c e unconditional is i s also equation f o r ADI stability of is the 37 scheme. 4 . 3 C o m p u t a t i o n a l D e t a i l s And We Sample P r o b l e m s have assumed t h a t t h e c o e f f i c i e n t s of as r e q u i r e d f o r the a n a l y s i s t o h o l d . stability various for e x t r a p o l a t i o n schemes f o r n o n l i n e a r w i t h known s o l u t i o n s and equations were s o l v e d u s i n g t h e ADI developed coding for this testing how to the ADI d i s c u s s i n g the evaluate o f t h e ADI considered, since c o n d i t i o n s defined at To see how The This i s centred intermediate equation The model Fortran we code u s e f u l when must In do time determine forcing functions addition, not the have boundary levels. Peaceman-Rachford formula i s c o n s i s t e n t w i t h t h e CN the into same as (6) coefficients scheme, of t h e ADI case, the functions difference accuracy coefficient d i s the 6 8 ) to evaluate The so a t t=n+* must ( 5 ) , w i t h a, b, c and about t = n + l / 2 . of equations. equations i s unchanged from the c o n s t a n t c o e f f i e n t s are The a l s o f o u n d t o be consistent. i t i s c o n v e n t i o n a l ( [ 9 ] , p. at t = n + l / 2 . which t. method. form to the scheme f o r b o u n d a r y v a l u e s our and series to incorporate v a r i a b l e c o e f f i c i e n t s scheme, c o n s i d e r of x, y and accuracy terms, a sample p r o b l e m s , remain constant, t o t e s t t h e e f f e c t of of s i m i l a r v a r i a b l e c o e f f i c i e n t s and equations requirement and was t h e p r o g r a m f o r t h e model Before ADI examine v a r i a b l e c o e f f i c i e n t s , and equations be To ( 5 ) are scheme assuming the known e x a c t l y . addition d i f f e r e n t problem. of a f o r c i n g f u n c t i o n t o an e q u a t i o n The forcing function itself must be poses a split in 38 some way t o m a i n t a i n the consistency of the ADI formula. C o n s i d e r i n g a s an example, u = u^ + u t we n o t e first (1-XD )u A + r (1-XD^)U and solving + F(x,y,t) , t h a t F ( x , y , t ) must be e v a l u a t e d a t t = n + l / 2 . 2 (9) y y * = (l+XD ^ + F ( x , y , (n+l/2)At)/2 2 = ( 1+XD. )u *+ F ( x , y , ( n + l / 2 ) A t ) / 2 w H 2 f o r u*** , we f i n d n+ that (10) u*** = ( 1 - X D ^ ) "MO+XD^u* + F where F A + , / 2- second s p l i t = Using F ( x ,y, (n +1 /2) A t ) . equation , W V V2) Substituting (10) i n t o the gives (l-XD. )(l-XD^)u 2 = O+XD^Ml+XD^Ju" + nH y ((l-XD ) + ( l + X D 2 E Y 2 < ))F^ V2 y Y = ( 1+XD^) ( 1+XD^)u* + F N + ^ w h i c h shows t h a t ( 9 ) i s c o n s i s t e n t . When solving time-dependent boundary t=n+*, a r e n e e d e d . to any p a r t i c u l a r averaging either time The v a l u e s u formula n + level could level, does n o t c o r r e s p o n d be approximated by a t t=n+1 a n d t = n , b u t t h e o v e r a l l a c c u r a c y o f t h e scheme w o u l d n o t be m a i n t a i n e d find equation, v a l u e s a t t h e i n t e r m e d i a t e time However, t h i s time. values Peaceman-Rachford * explicitly f o l l o w s from ([15]). I t i s possible to i n t e r m s o f t h e v a l u e s a t u " a n d u**' ; equation (7). Adding the split the equations gives ( ( 1 - XD - XD^) + (1 + XD^+ \L\ ) ) u« * = ( 1 - XD£- XD^) u ^ 2 + E Y y (l+XD^+XD^u' which i s j u s t u*** = ( ( 1 - X D ^ - X D e J ( )u n + , + (1+XD^+XEL )U )/2 . m x 1 1 + 39 We r e q u i r e u "* on nt the x - d i r e c t i o n boundaries, and i t is easily s e e n t h a t t h e d i f f e r e n c e o p e r a t o r s i n v o l v e v a l u e s o n l y on boundaries. This formula ensures that calculation b o u n d a r y v a l u e s p r e s e r v e s t h e o r d e r of a c c u r a c y We c a n now problems. ADI The has problem p r o b l e m s we the and coefficient a p p l y t h e ADI accuracy that, and ± predicted not as the levels, various absolute the above. error series of sample i n t e n d e d t o show t h a t accurate coefficient for variable t h e method a p p e a r s t o be the duy of Ax, at the conditions We chosen exact intermediate time are use calculated of the g r i d p o i n t s a t t i m e = 1.0 are given i n Table c = d = 1.0, solution u = exp(-x-y-t). according given to the 0((At) +At(Ax) ). 3 truncation 2 error in exact Table truncation We under error (I). of set The. the a = b = errors indicate, terms, whose h i g h e r d e r i v a t i v e s of e x p ( - x - y - t ) . which is 0.5 The vary method, which see t h a t t h e e r r o r s a r e s m a l l e r would c a n c e l l a t i o n of e r r o r 1 we for 1. problem, has as maximum For the c o n s t a n t c o e f f i c i e n t which to the consisting interior are interest. except results, v a r i o u s v a l u e s o f A t and r e s u l t s are at form the c o e f f i c i e n t s boundary The at + x problems boundaries where explained + cu y y s o l u t i o n u and s o l u t i o n on t h e the stable. + bu w of of t h e method. f o r the c o n s t a n t n o n l i n e a r problems, = au where t h e e x a c t and a sample p r o b l e m s a l l have the g e n e r a l (11) u produce to have chosen a r e although least conditionally The scheme those likely coefficients are than due is the to various 40 To test a variable coefficient a=0.5/(y+.1) , b=0.5/(x+.1) , 2 has under u=exp(-(x+.1)(y+.1)-t). ( I I ) and a r e much t h e same as For both these cases the e r r o r first ten checked for twenty f o r 200 time s t e p s ) , t h e t i m e and One treating the o t h e r as with a and c=-u</(u (y+.1)) approximated exact is coefficient of the Table 1 method As d e s c r i b e d i n C h a p t e r discretizing one The d=-Uy/(uy(x+.1)), same were table. derivative sample u , x as where is problem, u*, Uy are Uy a r e known e x a c t l y . f o r the v a r i a b l e 2, and equation f o r the v a r i a b l e c o e f f i c i e n t the the both model e q u a t i o n s c o n t a i n i s by e x t r a p o l a t i o n and solution stability i t s coefficient. and A by b as which r e s u l t s are given constant steps g i v e n i n the approximated d=l/(x+.1), becomes c o n s t a n t a f t e r indicating of the c o n v e c t i v e t e r m s . are solved space i s solved with s t e p s ( f o r some c a s e s , e r r o r s of the n o n l i n e a r terms the square these time The the case. to (11) c=1/(y+.1) and 2 true solution problem, The coefficient problem. Case I'll in approximating results u x with D y(1.5u*-0.5u " ). r t 1 0 accurate, by D 0 < u" contains and Uy approximations The the by results D yU and w 0 n A a l t h o u g h n o t as much b e t t e r is slightly as e x p e c t e d . later times, apparently caused by linking oscillating and by a p a r a s i t i c three time levels finally solution together. e x t r a p o l a t i o n d i d n o t b l o w up w o u l d For the and x 0 two-point e x t r a p o l a t i o n by c a s e IV D x (1 • 5 u - 0 . 5u ' ) runs w i t h t h i s e x t r a p o l a t i o n , however, the s o l u t i o n at found more several is unstable b l o w i n g up. This i s introduced numerically (That the one-point seem t o c o n f i r m t h i s . ) Since 41 the solution this p r o b l e m can space with At=.05, be AX=.10 d i d not a v o i d e d by a careful the two-point e x t r a p o l a t i o n two p r e v i o u s time l e v e l s , a s o l u t i o n We compared Taylor the series effect i n the test diffusion the errors we Here u = c a p p r o x i m a t e l y by The and solutions from be due to convergence made w i t h At solutions variable are the rate and for and the at Ax not (V). clear errors but only exact at less due and At. of the coefficient solution time=At is again was found examples except scheme i s obviously previous example, f o r the parasitic solutions. produce stable other to do not The the form r u n s were the reasons for of find These corresponding converge at At. To results. a c c u r a t e than the c h a n g e s o n l y as seem t o be variable t h i s problem, several c o e f f i c i e n t p r o b l e m , but e x p e c t ; the the found 2 A l l the As p r e s e n c e of chosen to Ax a a=.5u/((y+.1) u), l a t e r t i m e s ; the stable. for and by expansion. g i v e n under b l o w up solution setting d as necessary. coefficients solution series is t e s t v a l u e s of (12) The longer unconditionally t h i s may not time this solution exact nonlinear n Taylor AX=.10 of 1 . Su^-0. 5u "' r e s u l t s are At=.05, f o r the and u=exp(-(x+.1)(y+.1)-t). time=At approximating solve and 2 problem. of effect terms, b=.5u/((x+.1) u) the of requires for expansion w i t h u s i n g the difference To no choice i t appears steps. Since no b l o w up, the rate this we are nonlinear problem. We and the conclude from t h e s e sample problems t h a t methods chosen for handling nonlinear the ADI terms scheme produce 42 acceptable stability solutions. The larger errors and f o r t h e n o n l i n e a r e q u a t i o n s a r e due specific solution and the fact that in the conditional part to the equations have c o e f f i c i e n t s which c o n t a i n exact s o l u t i o n s or exact derivatives. The these will error that results always necessarily be noted be in the case as from t h e a p p r o x i m a t i o n s to t h e same d i r e c t i o n , w h i c h i s c e r t a i n l y not f o r t r u l y nonlinear equations. w e l l t h a t the e r r o r produced model equations are highly d i f f u s i v e ; will have o n l y a s m a l l c o n t r i b u t i o n It should by e x t r a p o l a t i n g t h e c o e f f i c i e n t s of the c o n v e c t i v e terms i s of l i t t l e the terms concern since the c o n v e c t i v e terms t o the o v e r a l l solution. At Ax .05 .05 1.5E-5 4.1E-6 5.5E-4 3.3E-4 1 .6E-4 .05 .1 0 5.8E-5 1.5E-5 5.7E-4 3.4E-4 1.3E-4 .1 0 .05 1.9E-5 5.9E-6 1.1E-3 1.3E-3 2.1E-3 .10 .10 6.2E-5 1.6E-5 1.1E-3 6.4E-4 2.9E-4 .20 .20 2.5E-4 6.7E-5 2.3E-3 4.8E-3 3.9E-4 Table 1. (I) (II) final (9) a r e (9) (V) problems. Equations f o r m o f t h e ADI scheme f o r t h e m o d e l e q u a t i o n s i s f o u n d by a p p l y i n g t h e p r o c e d u r e s equations (IV) R e s u l t s o f ADI s o l u t i o n s o f s a m p l e 4.4 A p p l y i n g ADI To The M o d e l The (III) and discussed in (10) o f C h a p t e r 2. this The s p l i t chapter to equations f o r 43 (l-Xa D - X a D ~ JH"** 2 2 (13) 4 y ^ [ 3 and f o r (10) a r e O-XbzD, -Xb„D (14) ^ 2 0y 3 where t h e c o e f f i c i e n t s , e>t )S*+* )S* f, a t time=n+1/2. = (l+Xa D 2 +Xa ^ )H + Atf,/2 n 2 0 y + Xa D , 4 < ) H + * A t f , /2 n1 y = ( l + X b ^ + Xb D . ) S* + A t f / 2 * 2 3 y ( 1 - X b ^ -Xb D 2 ^ ( l - X a ^ - X a D ^ )U** and e v a l u a t e d = ( 1 +Xa , D M 01 2 = ( 1 + X b D ^ + Xb D^ J S * ^ A t f / 2 and f 2 2 1) a r e d e f i n e d as i n 2 Chapter 2 44 V. In the numerical first part t e s t s of the involve conserving. sets was next of We these i n i t i a l of input be is and A the model, both w i t h improvements t o the ADI equations and as m e a s u r e d by using The with time series the plotted of continued taken to more m o d e l l e d and with methods u s e d . nonlinear approximations mass, as m e a s u r e d by H, product the of H and specified s e t t i n g the velocities x-direction boundaries salinity initial values u and and to the total We but the model salinity, consider impose of H and S v to zero the closed on on all the r e s p e c t i v e l y , i t i s c l e a r that i n the conditions S. domain, b o u n d a r i e s , and For of w h i c h , i t various p h y s i c a l system being constant total both manner. w h i c h m i g h t be maintaining mass and of recommendations f o r By and are these t e s t s i s to ensure t h a t a p p l y i n g the the over boundaries. closed approximate s o l u t i o n s interpretation steps numerical of conserve equations to produce a Tests purpose scheme the They t e s t i n g the model e q u a t i o n s .physical work on The to conclude with represent discussed. straightforward boundary c o n d i t i o n s We are preliminary computational d e t a i l s in a devoted follows. Preliminary the several boundary c o n d i t i o n s , n e i t h e r results 5.1 the applied values. accurately equations to v e r i f y that and found, could t h i s chapter, boundary c o n d i t i o n s also discuss initial section of model s i m p l i f y i n g the s y s t e m , a l l o w i n g us CHAPTER 4 we system should use constant remain values ythe constant. of H and S 45 everywhere is except at a single grid point. At t h i s p o i n t , i n the c e n t e r of the domain, H or S or both larger than constant the constant Three i s set to a value runs were v a l u e s o f 20 m e t e r s f o r H and 2 p . p . t . s i z e s o f Ax=Ay=At=0.1. Mass field. and HS are The r e s u l t s a r e g i v e n calculated simply which made with f o r S, a n d s t e p i n Table 2 below. by summing o v e r the g r i d points. H(6,6) t =0 S(6,6) t=0 Mass t =0 HS t=0 2.0 2422.0 4844.0 2420.2 4840.4 22.0 2.5 2422.0 4855.0 2422 . 9 4840.7 28.0 2.5 2428.0 4870.0 2421 . 4 4840.9 2. R e s u l t s of C o n s e r v a t i o n The method a p p e a r s t o c o n s e r v e an the f i r s t crude error over 100 t i m e c a s e , and s l i g h t l y Tests both q u a n t i t i e s s t e p s of l e s s than larger the derivatives point i n terms of c o n s e r v a t i o n . well, .1 p e r c e n t f o r For the The r e s u l t s o f t h e i n d i c a t e , not s u r p r i s i n g l y , that the sources, are approximated, very f o r t h e o t h e r two. 11 by 11 g r i d u s e d , t h i s i s v e r y g o o d . s e c o n d and t h i r d t e s t s a l s o larger HS t=1 00 22.0 Table with Mass t=1 00 the less accurate and t h e l a r g e r the the overall spatial error 46 5.2 Initial C o n d i t i o n s At The For purposes a s s i g n e d as initial layer, and H, of time Here test discharge 2 ppt f o r S, were the point problem larger than the initial f u n c t i o n at the solutions The iterations, and f o r H and S would problem equation The y stems = -(fHS y near-discontinuity, from S d e p e n d e d on d o m a i n , and the a value nearby, we the of we S at the used a weekly resulted in of extremely f o r b o t h H and S at was no indication that the having to find the s o l u t i o n of 1, w h i c h i s + 2kSH x + kHS^)/2fS . c r e a t e d by a l a r g e i n p u t To a v o i d the effect value, of this assume [ 7 ] t h a t S a t t h e d i s c h a r g e point arbitrarily along discharge These r e s u l t s p e r s i s t e d d e s p i t e of l a r g e c h a n g e s i n s a l i n i t y direction We river i m m e d i a t e e f f e c t on p o i n t s i n t h e r e g i o n ; We upper converge. (1) p r o d u c e s l a r g e v a l u e s of H. region. the v a l u e s c h o s e n were at the level there t o t h e l a r g e g r a d i e n t s of S impact were nearly-discontinuous nature time (14), from Chapter (1) H an by those g r i d p o i n t s near the d i s c h a r g e . has The v a l u e s of the s o l u t i o n s a t time At several values the h e i g h t of a s s i g n e d t o H and created to f i v e . Due except field testing. s e r i e s with S equal unlikely constant t h e s a l i n i t y d e f e c t , S. t h e v a l u e s we s e r i e s we To testing, c o n d i t i o n s to both 20 m e t e r s f o r H and point. Land Boundary chose land varied S linearly defect five boundary, over is, the i s smoothed o v e r the grid two that points in either grid points into the these p o i n t s . a l s o want t o f i n d a v a l u e of H a t the discharge point 47 such that H S. To f i n d 2 is H S is assumed point, for the discharge includes immediately U in e q u i l i b r i u m with nearby values s u c h an H , we use that but points initially 2 the the following region region procedure consists up t o north and-south. but of not only [7], of It one grid i n c l u d i n g the grid The l a n d b o u n d a r y c o n d i t i o n is U/C = - ( f ( H S ) + k(H S) ) 2 with U given gradient Coriolis at of the H S force, x discharge across 2 point the and region and we a p p r o x i m a t e it (H SL. = (H SL - 6(H S) (H S) = (H S)< + 5(H S) 2 , 2 y 2 U = 0 is elsewhere. negative, The due t o the using 2 (2) 2 where the K, the the south K-1. discharge first 2 y grid find 5(H S) that a = 5(H S)/Ax 8(H S), we 2 K , and denoted by t h e n o r t h by K+1, and t h e increase second subscript first i n H S from p o i n t 2 order to K to approximation to 2 we discharge that U = 0 at all grid points except 2 Denoting can point further point x the approximate 2 discharge note = k(H S) -f((H S)<*> we assume . so these p o i n t s . (H S)©« , If the the 2 y 2 the is 2 -f(H S)y at to is is 2 point 2 (LAx,KAx) point from (2) (H S) To point by K - 1 . We see (H S) 2 K r J i n i t i a l constant -f(H S) 2 on the field value by b o u n d a r y n o r t h of by - (H S)^, 2 that and a t (H S) the )/5Ax X has five grid 2 . the points same value at n o r t h and s o u t h , the we 48 have k(H S) 2 and = -f((H S)eo so a t t h e d i s c h a r g e (H S) 2 )/5Ax ( C + | point, = f6(H S)/Ax -U/C " 2 x - f((H S)co 2 - (H S) )/5Ax 2 2 | C + ) or (3) - U A x / ( C f ) = 1 . 2 8 ( H S ) - 0 . 2 ( ( H S ) o o 2 We a l s o r e q u i r e t h a t adjacent point in t o the discharge i s ((L-1)Ax,KAx), the x-direction " 2 the grid point p o i n t be e q u a l (H S) ) . 2 K inside t h e domain t o (H S)o© . Since 2 H S c a n be e x p a n d e d i na Taylor 2 this series about LAx, which i s , t o O(Ax), H S((L-1)Ax,KAx) = H S(LAx,KAx) - Ax(H S(LAx,KAx)) 2 2 2 or (H S)oo (4) = (H SV 2 2 We know S e v e r y w h e r e a l o n g and ( H S ) point, in equations with are t h e boundary, , our i n i t i a l o e " 2 (H S) 2 U at I C - 6(H S)) 2 2 K the discharge equations s o l v e d , a n d once 6 ( H S ) a n d ( H S ) ^ 2 H(LAx,(K±i)Ax), i = 1 , . . . , 5 linear first from (H S)(LAx,(K±1)Ax) 2 The values result of a l l this 2 maintained. region so that sufficiently 2 Finally, determined using a from ( 2 ) . t o smooth l a r g e of finding H initial at the an e q u i l i b r i u m i n ( H S ) When i n c o r p o r a t e d i n t o t h e m o d e l , smoothed t h e r e s u l t s by a s s u m i n g H S 2 method initially determined f o r (H S) , the values f o r i s simply a . 2 i=1,...,5, a r e e a s i l y of S, a n d t o p r o v i d e discharge , a r e found to (H S) order d i f f e r e n c e approximation H((L-1)Ax,(K±i)Ax), unknowns Two known, H ( L A x , ( K + 1 ) A x ) a n d H ( L A x , ( K - 1 ) A x ) c a n be (1). . 2 c o n d i t i o n , so t h e o n l y (3) and (4) a r e 6(H S) and ( H S ) . two unknowns a r e e a s i l y from is 2 + (f/5k) ( (H S)o* 2 this t o produce reasonable is procedure results. 49 5.3 F i n d i n g H On In The 1, Chapter conditions Boundaries we mentioned f o r H were p r e s c r i b e d on boundaries. The computational only both that the details the boundary ocean are and given land in this section. On the reasonable and open to t h a t any (and H on two S) time current modified by on line. a grid at point For the H at the gradients Except at used using the H on the initially, e x t r a p o l a t i n g from time l e v e l , the river The l a n d b o u n d a r y , we direction. discharge the When boundary, the an point. land boundary, solving s i g n s of p r o c e s s of of d i f f i c u l t i e s the of first we values two nearest time at the formula f o r H has boundary H calculate extrapolation the is is been further interior points use f o r example. s e c t i o n w i t h as m i d d l e of number to It change r a p i d l y , values on the ocean more c l o s e l y what i s o c c u r r i n g w i t h i n t h e H on In the unknown. A f t e r the model e q u a t i o n T h a t i s , we (iAx,0), previous the values e x t r a p o l a t i n g from the thus r e f l e c t is due domain. levels i n C h a p t e r 2. at the occur are the i s set to the previous solved within o c e a n b o u n d a r i e s by discussed H assume, h o w e v e r , t h a t H d o e s n o t where H the boundaries c h a n g e s t h a t do hence level, ocean it must s o l v e e q u a t i o n initial Since (l) for value the the be domain. (1) from input value the d i s c h a r g e must boundary point solved is in in each s o u t h e r n p o r t i o n of y - d e r i v a t i v e s must be of the reversed. finding a solution to equation were e n c o u n t e r e d . One was the (1), a e f f e c t of 50 the salinity discharge S that defect point. were values large Even w i t h section, this much l e d to a large gradient, points a point. Integration evaluating a input Ax, cubic spline values to i n the l a s t to solve near the Ax a s a s t e p size results. t o u s e an i n t e g r a t i o n which river Sy, i n e a c h d i r e c t i o n a n d methods u s i n g requires between our g r i d p o i n t s . fitting zero than the boundary and i n t e r i o r t h a t was d i f f i c u l t i s necessary therefore smaller at the smoothing procedure d e s c r i b e d produced p h y s i c a l l y impossible It input compared t o a d j a c e n t produced a s t i f f equation only as F o r t e s t p u r p o s e s , we a s s i g n e d points. discharge chosen values F o r S, t h i s through S size f o r S and U a t i s accomplished at the g r i d the curve at the required p o i n t s . at the g r i d points corresponding step For p o i n t s and U, which point (JAx,(K+1)Ax) and ( J A x , ( K - 1 ) A x ) t o t h e s o u t h , the north n e e d some f u n c t i o n w h i c h d e s c r i b e s exponential curve with U(JAx,KAx) i s the input U(JAx,(K-1)Ax) are a (JAx,KAx) negative discharge very nearly and i t s behavior. large We the nonzero values of U i n equation t e r m - u / ( 2 f S C H ) on t h e right side points chose exponent, such value zero. grid and This, an that i t i s hoped, i s a ( 1 ) , we must of we U(JAx,(K+1)Ax), r e a s o n a b l e a p p r o x i m a t i o n t o U a t t h e " r i v e r mouth". these is t o the l a n d boundary, but nonzero between t h e d i s c h a r g e to by To include reintroduce ( 1 ) , and the ODE becomes (5) Hy = -(U/(CH) + f H S A terms. further y +2kSH x + kHS )/(2fS) x p r o b l e m a r i s e s when e v a l u a t i n g Approximating H and S with the x - d e r i v a t i v e the second order difference 51 equation detailed earlier the domain. same Unfortunately, time l e v e l values must (say estimated previous time interior p o i n t s , we The boundary from the have be two interior a p p r o x i m a t e d , and boundary v a l u e s subroutine, f o r H and equation to which differentiation corrector. H on an formula The land solved. the (5) H This for; boundary at the S. extrapolated and S we can now in new results using the [10] from f o r H and extrapolated method consult the Once s o l v e d , i s found u s i n g iterative should at (5) values u s u a l way. as a p r e d i c t o r and reader the S a t t i m e n+1, (5) equation is on a l l boundaries are for H which replace solution s o l v i n g equation approximate l e v e l s i n the values inside p o i n t s must be s o l v e the model e q u a t i o n s for previous points Rather than e x t r a p o l a t i n g provide first values interior A to are p o i n t s or i n some way. levels at g r i d these i n t e r i o r t i m e = n + l ) we e i t h e r at the be involves values values. t h e UBC a GEARB backward c h o r d method as a f o r a more d e t a i l e d explanat ion. We can n o t e h e r e t h a t s i n c e H must be change r a p i d l y making the e x t r a p o l a t i o n wish to avoid intermediate time. for e x t r a p o l a t e d , and We H first having to find land time l e v e l , which does not do this by splitting the boundary are thus only equation required accurate, values c o r r e s p o n d t o an i n the y - d i r e c t i o n , then i n the boundary v a l u e s less since H we at the actual s u c h t h a t we solve x-direction. Land f o r t i m e s n and n+1. 52 5.4 Boundary C o n d i t i o n For The no-flow to boundary Using d e r i v a t i v e has matrix condition c o n d i t i o n , which zero. a one equations salinity defect particularly first second condition is the first For As rapidly the on we land boundary, such c a s e s , it is very on the important to with at least t h e same accurate usually approximation centred about the this formula to grid i s a reasonable the no-flow point on domain. approach, no p h y s i c a l m e a n i n g . involving only i n t e r i o r On b u t on Instead points; the the the we f o r the boundary, f o r example, t h i s i s (6) S.j - 4S j /3 + S^ /3 t> = 0 , j = e q u a t i o n s a p p l y on t h e o t h e r Incorporating t h a t e a c h row elements system (6) corresponding to the r i g h t is lost. performed to 1,...,L. boundaries. i n t o t h e ADI matrix equations the boundary Rather than try of to contains two, s o l v e a system a on e a c h o f t h e s e rows b e f o r e s o l v i n g example, the f o r S, of the d i a g o n a l ; the t r i d i a g o n a l i t y h a l f - b a n d w i d t h of t h r e e i n s t e a d For have n o t e d , d i s c h a r g e r e g i o n (depending l a n d b o u n d a r y , s u c h a v a l u e has see the scheme. open o c e a n b o u n d a r i e s , Similar to the tridiagonal. boundary, r e q u i r i n g a value f o r S o u t s i d e the western derivative remain order use a d i f f e r e n c e d e f e c t i s the the boundary c o n d i t i o n here as t h e ADI A salinity m a j o r a d v a n t a g e ; when i t i s d i s c r e t i z e d , input values used). accuracy the order accurate approximation the r i v e r the for requires setting changes near approximate S first two row of rows of t h e m a t r i x e q u a t i o n two the with a reduction the whole we is system. for S at 53 t=n+* a r e 1 -4/3 1/3 a1 a 2 2 0 0 3 T h i s c a n be r e d u c e d t o (1-a ,/(3a )) A the 22 2 1 Reducing each such returns (-4/3-a /(3a )) 2 3 2 3 row system f o r both 0 23 2 2 A left 2 3 and r i g h t to tridiagonality, and boundaries increases the o p e r a t i o n c o u n t by o n l y 0 ( L ) . 5.5 T e s t s Of The M o d e l Incorporating conditions Equations these methods equations. we c a n p r o c e e d w i t h t e s t s Input to the midpoint point are input. e q u i l i b r i u m with these Input equilibrium i n section values S=5.0 p p t , value F o r each time 4.3 u, so experimental u=0.3 m/sec, and t h e i n t o t h e model a t time t h e f i n a l time at the located step, the corresponding as t o be i n run c o n s i s t s of a pulse step, t=50. corresponding Initial a r e S=2.0 p p t a n d H=20.0 m t h r o u g h o u t t h e d o m a i n . same u n t i l model of S and u a t time t=0. f o r H, w h i c h i s H=17.66. values are introduced point section The t h i c k n e s s o f t h e u p p e r l a y e r , H, f o r the f i r s t with the f u l l t o be a s i n g l e g r i d o f t h e x=LAx b o u n d a r y . i s c a l c u l a t e d as e x p l a i n e d inflow, of d e f e c t , S, a n d x - d i r e c t i o n v e l o c i t y , this boundary t o the program c o n s i s t s o f t e s t data i n f l o w , w h i c h was c o n s i d e r e d salinity and i n t o t h e FORTRAN p r o g r a m t e s t e d i n t h e f i r s t of t h i s c h a p t e r , at for initial conditions The inflow t=0, and remain t h e 54 Based on experimental this the results of the first run, run i s i d e n t i c a l w i t h the first until the second t=l5. At p o i n t the boundary c o n d i t i o n s a t the i n f l o w a r e changed t o u=0 a n d dS/dx=0, c o r r e s p o n d i n g fresh water into t h e two p r e v i o u s continued until to shutting the model. time time steps H i s found in the usual off the supply of by e x t r a p o l a t i n g f r o m way. This run is s t e p 40. 5.6 I n t e r p r e t a t i o n Of R e s u l t s The r e s u l t s a r e p r e s e n t e d contour the i n F i g u r e s 3 t o 35 a n d c o n s i s t o f p l o t s of H and S and a c o n t o u r pressure p l o t of H S 2 (representing f i e l d ) o v e r l a y e d w i t h a v e c t o r p l o t of v e l o c i t y . The v e c t o r p l o t i s normalized to the time level, whose v a l u e i n m/sec point source of f r e s h water i s a t y - a x i s , a t x=600 km, y=300 largest velocity at the i s p r i n t e d on t h e f i g u r e . the midpoint of the The right km. Due t o t h e a l g o r i t h m u s e d f o r c o n t o u r i n g , a number o f p l o t s have r e g i o n s of unusual o c c u r s have v a l u e s contouring cannot produce f o r the western ignored; the accurate contours The by used, these around the i n f l o w p o i n t . over the next few t i m e where t h i s t h e same, a n d accurate contours. half o f t h e domain the On a l l should be a r e a l w a y s smooth c u r v e s . F i g u r e s f o r t=2 r e f l e c t the i n i t i a l large gradients procedure The a r e a s that are a l l approximately program p l o t s the r e s u l t s looking contours. of gradients numerical d i f f i c u l t i e s S. Despite induce the created smoothing unlikely v a l u e s of H By u s i n g a l a r g e number o f iterations steps, values f o r H a r e smoothed, as can 55 be s e e n by The r e g i o n of low Coriolis system. the t=5. f o r c e , growing This along in this f o r the p r e s s u r e be p r i m a r i l y contours 2 Figures layer time step for later t h i c k n e s s and along the s u r f a c e . term HS 2 reflect By time conditions "reflected" i n t o the domain. for this H, w h i c h may show t h e not The flow to l i e exactly is the flow w e l l e s t a b l i s h e d and i n c r e a s e i n the by a d v e c t i o n of t h e upper f r e s h water i n [ 6 ] , i n c r e a s i n g a m o u n t s of f r e s h s o u t h a s w e l l due to d i f f u s i v e l i m i t s of t h e d o m a i n a n d flow. of t o have an a f f e c t I t appears that the effects. the on flow the n o r t h e a s t c o r n e r , c a u s i n g A p o s s i b l e cause i s the g i v e too l a r g e a value boundary both upper is being f l o w back extrapolation f o r H on scheme the boundary i n area. While the results c o n d i t i o n s were u s e d u n t i l increases along change v e r y Being for pattern a spreading from to contours. i n force there begin away velocity s t e p s show o n l y an s t e p 30, t h e l a y e r t h i c k n e s s and and t h e e f f e c t s of f r i c t i o n ; this As n o t e d w a t e r move west and due i n c r e a s e i n the t h i c k n e s s T h a t t h e v e l o c i t i e s do 15 time north, r e g i o n ; to the south, H d e c r e a s e s . i s p u s h e d t o t h e l e f t of t h e s e By slowly l a r g e r as more f r e s h w a t e r f l o w s i n t h e northwards. the H S moves i s a c c o m p a n i e d by an upper l a y e r figures salinity are not time s t e p 100. the n o r t h e r n little boundary, but from those identical to the the second e x p e r i m e n t a l shown, a t time first step the The same boundary reflecting otherwise flow the results t = l 5 , the results 40. test u n t i l run begin a t time step 20 (Figures 56 27, 30 and "turned their 33). A t t=16, off"; expect the v a l u e s t o s l o w l y r e t u r n to i n i t i a l conditions. As c a n be is in due t o a d v e c t i o n and into to 3.6 the next ppt the until entrainment while that the of H to H layer f r o m one of HS 2 the time S is discussed flow continues inflow point, to about experimental runs the HS left of just the show t h a t t h e n u m e r i c a l model sign domain. of any instabilies. r u n a l s o shows t h a t a g i v e n p h y s i c a l s i t u a t i o n , other first with solving solutions, so no same, b u t time S was v a l u e s o f H were t h o s e models the e f f e c t on to illustrate Both available. the first test T h e r e a f t e r , the v e l o c i t y slowly decreased found to t h e model e q u a t i o n s . f i g u r e s are s t e p 40. the modelled. used the boundary c o n d i t i o n s of above u n t i l remained the i s adequatedly t e s t s were c a r r i e d o u t w h i c h s e r v e difficulties 2 to r e s u l t s w i t h no due by the ( F i g u r e 35) perhaps i n the n o r t h e a s t c o r n e r of the produced unstable The The the of system r e t u r n i n g to r e s t , further about d e f e c t moves s o u t h p u s h e d south produces reasonable Two step m. T h e r e a p p e a r s t o be an eddy second water f r o m a h i g h of contours. these this decrease f l o w , whose s p e e d s ( t h e maximum of w h i c h i s g i v e n on Both The lower value decreased r e g i o n of h i g h e s t s a l i n i t y "reflection" The maximum has figures) diminish rapidly. the of s a l t i e r These v a l u e s decrease a t t=40 t o a b o u t 22 The s e e n i n F i g u r e s 27-35 what h a p p e n s , a s t h e g r a d i e n t s of S and the upper l a y e r . 28 m. f r e s h w a t e r f l o w i n t o t h e model i s would fact we the be in over time to equilibrium. t h e s y s t e m o f a f l o w of v a r y i n g u S=2. This salinity. 57 The as 0.1 salinity ppt per developed a l a r g e number of t i m e s other conditions, unsuccessful but with i n d i c a t e d above. time steps. diverging a the model little quickly ( s o o n f o l l o w e d by showed t h e H). solution for S to most quarter within two of or i t e r a t i o n s again that three led to a seem grid The the to be r e s u l t i n g values any section. caused first by and should numerical be made t h e two tests first. In by both cases, the l a r g e g r a d i e n t s near , i t r e s u l t e d from reducing ppt after input s a l i n i t y i n equation gradients ( 5 ) , the the c o a s t e n t e r i n g the model at the reaching defect dropped producing o f S on The c h a n g e s can be d i s c u s s e d last f r o m 5.0 points instabilities will of defect state solution. surrounding the model. problem i s i l l u s t r a t e d In the equations have c h o s e n t o s o l v e them t h a t when i m p r o v i n g so t h e y a t t h e end salinity t h a n any boundary Proceeding obvious c o a s t a l boundary. created resulted must be d e a l t w i t h b e f o r e instabilities input usual s p a t i a l mesh s i z e one methods we the e q u a t i o n s , discussed the s e v e r a l s h o r t c o m i n g s i n the model be c o n s i d e r e d difficulties The involved solution. numerical initially test Instabilities Recommendations For the the of S I n c r e a s i n g t h e number of There are in i n the v a l u e s rate r a t e of as actually diverging. The 5.7 decreased at the t i m e s t e p ; even a t t h i s instabilities Iterating be d e f e c t was land a the near-steady below that S at apparently boundary increase to values inflow. the ODE; larger 58 Results results of from the from a decrease second indicate i n mesh s i z e . t h a t t h e same Since the right-hand side ( 5 ) d e p e n d s on x - d i r e c t i o n d e r i v a t i v e s o f H a n d S approximated which by d i f f e r e n c e e q u a t i o n s , t h e s e v a l u e s become (assuming a d j a c e n t v a l u e s o f H and S a r e a p p r o x i m a t e l y and gets (5) even s t i f f e r as the r e s o l u t i o n t h e s o l u t i o n b l o w s up a l o n g increasing the fourth order. apparent, the accuracy the of coastal the A method o f s o l v i n g t h i s p r o b l e m numerical problem segments of as when readily scheme u s e d f o r IVP, i s c r e a t e d when boundary and to now the contains south v a l u e problem of i t seems l i k e l y model, a realistic involving number of solving present further processes some conditions for H and s o r t of h i s t o r i c a l also overlay a geostrophic current f i e l d , model i s p r e s e n t l y c o n s i d e r e d a t r e s t ) , S both the First, are average. ends. of a d d i t i o n s s h o u l d be c o n s i d e r e d . s e t of i n i t i a l perhaps the p h y s i c a l ends in p r o b l e m s a s a BVP, p r o b a b l y h a v i n g b o u n d a r y l a y e r s a t b o t h To more a c c u r a t e l y r e f l e c t the ( I V P ) becomes encountered that i tw i l l To Between (BVP), as t h e v a l u e s f o r H a t both Considering the d i f f i c u l t i e s an not ( 5 ) c a n be s o l v e d a s u s u a l . segments, however, t h e i n i t i a l a boundary v a l u e problem (5) is as y e t u n r e s o l v e d segment segment, e q u a t i o n are f i x e d . even t h e c o a s t where r u n - o f f e n t e r s t h e m o d e l . the n o r t h of the northern these Again d e r i v a t i v e approximations to the t r u e domain i s c o n s i d e r e d ; t h e l a n d southern t h e same) c o n d i t i o n s m i g h t be c o n s i d e r e d . Another two boundary, are larger increases. a l t h o u g h a v a r i a t i o n of t h e a v e r a g i n g initial problem needed, We initially and a t each time a should (the level. 59 As n o t e d i n the previous section, reflecting particularly behavior on some method o f more a c c u r a t e l y the open-ocean boundaries a t t h e i n t e r s e c t i o n of t h e c o a s t a l boundary. the model e q u a t i o n s s h o u l d i n c l u d e a d v e c t i o n c u r r e n t wind f o r c i n g , and w i n d - i n d u c e d entrainment i s required, term. Last, terms and m i x i n g s h o u l d be a d d e d t o W , t h e The e f f e c t s o f t h e s e c h a n g e s on t h e n u m e r i c a l methods f o r t h e s i m p l i f i e d model w o u l d have t o be d e t e r m i n e d . 60 Figure 4 - Test 1 : S at time step 5 61 0.0 Figure 100.C 200.0 5 - Test 300.0 KM . 1 : S at 400.0 soo.o soo.o time step 10 time s t e p 15 KM Figure 6 - Test 1 : S at 62 KM -200.0 Figure 8 - Test 1: S at too. a soo.o time step 30 63 Figure 10 - T e s t 1: S at time step 50 64 ^P- i 0.3 Figure t 0.0 Figure 3 I 100.3 200.3 300.3 1 - 1 200.3 300.3 KM 11 - T e s t I 100..3 ! 230.0 12 - T e s t 400.3 1: H a t : i 300.0 KM 1: H a t SOO.O SOof'3 I i ; 400.0 SOO.O SOO.O time 1 400.0 time 0 step 2 1 SOO.O step 5 ~ 0 S00.3 65 Figure 14 - T e s t 1: H a t time s t e p 15 66 Figure 15 - T e s t Figure 16 - Test 1: H a t time step 20 H at time step 30 67 Figure 18 - Test 1: H at time step 50 68 Figure 20 - Test 1: H S 2 at time step 5 69 Figure 22 - Test 1: H S at 2 time step 15 70 71 Figure 26 - Test 1: H S at 2 time step 50 72 F i g u r e ' 28 - Test 2: S'at time step 30 73 I 3.3 Figure I lttD.3 29 - i I 2M.3 Test 330.3 KM 2: S at l 40D.3 time i Sflfl.3 step o V 503.3 40 74 n 0.0 Figure i 100.0 ! 200.0 30 - T e s t 1 300.0 KM 2: H a t 1 1 r 400.0 SOO.O time step 20 time s t e p 30 KM Figure 3 1 - T e s t 2: H a t o 500.0 Figure 32 - Test 2: H a t time s t e p 40 76 Figure 34 - Test 2: H S at 2 time step 30 77 sac.s 9.3 •200.3 300.3 300.3 saol2 / 600.3 KM Figure 35 - T e s t 2: H S at 2 time step 40 78 BIBLIOGRAPHY 1. Bo P e d e r s e n , F., A M o n o g r a p h on T u r b u l e n t E n t r a i n m e n t a n d F r i c t i o n i n Two L a y e r S t r a t i f i e d F l o w , I n s t . Hydrodynamics and H y d r a u l i c E n g e n e e r i n g , Tech. Univ. Denmark, L y n g b y . S e r i e s P a p e r No. 25. (1980) 2. F o r s y t h e , G.E. a n d Wasow, W.R., F i n i t e D i f f e r e n c e M e t h o d s f o r P a r t i a l D i f f e r e n t i a l E q u a t i o n s , J o h n W i l e y & S o n s , New York ( 1 9 6 0 ) . 3. G e a r , C.W., N u m e r i c a l I n i t i a l V a l u e P r o b l e m s i n O r d i n a r y D i f f e r e n t i a l E q u a t i o n s , P r e n t i c e - H a l l , Englewood C l i f f s , W\~J~. ( 1 971 ) . 4. G o t t l i e b , D. and G u s t a f s s o n , B., G e n e r a l i z e d Du F o r t - F r a n k e l Methods f o r P a r a b o l i c I n i t i a l - B o u n d a r y V a l u e P r o b l e m s , SIAM J . Numer. A n a l . , 13 (1976) 129-144. 5. K r i e s s , H. and O l i g e r , J . , M e t h o d s o f t h e A p p r o x i m a t e S o l u t i o n o f Time D e p e n d e n t P r o b l e m s , GARP P u b l i c a t i o n s S e r i e s No. 10 (1973) . 6. L e B l o n d , P.H., Emery, W.J., N i c o l , T.O., A M o d e l o f R u n o f f Driven C o a s t a l C i r c u l a t i o n , Submitted f o r p u b l i c a t i o n t o C o a s t a l , E s t u a r i n e a n d S h e l f S c i e n c e , (1983) 7. L e B l o n d , P.H., D e p a r t m e n t o f O c e a n o g r a p h y , U n i v e r s i t y o f B r i t i s h C o l u m b i a , V a n c o u v e r , B.C., p e r s o n a l c o m m u n i c a t i o n . 8. L e e s , M., An E x t r a p o l a t e d C r a n k - N i c o l s o n D i f f e r e n c e Scheme for Q u a s i l i n e a r P a r a b o l i c E q u a t i o n s , Numerical S o l u t i o n of P a r t i a l D i f f e r e n t i a l E q u a t i o n s , Ames, W.F. (Ed.), Barnes & N o b l e , New York ( 1 9 6 9 ) , 193-201. 9. M i t c h e l l , A.R. a n d G r i f f i t h s , D.F., The F i n i t e D i f f e r e n c e Method i n P a r t i a l D i f f e r e n t i a l E q u a t i o n s , John W i l e y & S o n s , C h i c h e s t e r (1980) 10. M o o r e , C.F., I n t e g r a t i o n o f S t i f f S y s t e m s o f O r d i n a r y D i f f e r e n t i a l E q u a t i o n s H a v i n g a Banded J a c o b i a n , C o m p u t i n g Centre, U n i v e r s i t y of B r i t i s h Columbia (1979)., 11. M o r r i s , J . L . , On t h e N u m e r i c a l S o l u t i o n o f a H e a t E q u a t i o n A s s o c i a t e d w i t h a T h e r m a l P r i n t - H e a d , J . Comp. P h y s . , 5 ( 1 9 7 0 ) , 208-228. 12. R i c h t m y e r , R.D. a n d M o r t o n , K.W., D i f f e r e n c e M e t h o d s f o r I n i t i a l - V a l u e P r o b l e m s , I n t e r s c i e n c e P u b l i s h e r s , New Y o r k (1967) . 79 13. R o y e r , T . C , On t h e E f f e c t o f P r e c i p i t a t i o n and R u n o f f on C o a s t a l C i r c u l a t i o n i n the G u l f of A l a s k a , J . Phys. Oceanogr., 9 ( 1 9 7 9 ) , 555-563. 14. R o y e r , T.C., H a n s e n , D.V. and P a s l i n s k i , D.J., C o a s t a l F l o w i n t h e N o r t h e r n G u l f o f A l a s k a a s O b s e r v e d by Dynamic Topography and S a t e l l i t e - T r a c k e d Drogued D r i f t Buoys, J . Phys. Oceanogr., 9 (1979), 785-801. 15. R o y e r , T . C , C o a s t a l F r e s h Water D i s c h a r g e i n t h e N o r t h e a s t P a c i f i c , J . Geophys. R e s . , 8 3 , C3 (1982) 2017-2021. 16. S o m m e i j e r , B.P., Van Der Houwen, P . J . a n d V e r w e r , J.G., On t h e T r e a t m e n t o f T i m e - D e p e n d e n t B o u n d a r y C o n d i t i o n s i n S p l i t t i n g Methods f o r P a r a b o l i c D i f f e r e n t i a l E q u a t i o n s , Int. J . Numer. M e t h . E n g . 17 (1981> 335-346. 17. V a r a h , J.M., S t a b i l i t y o f D i f f e r e n c e A p p r o x i m a t i o n s t o t h e Mixed I n i t i a l Boundary V a l u e Problems f o r P a r a b o l i c S y s t e m s , SIAM J . Numer. A n a l . , 8 ( 1 9 7 1 ) 598-615. 18. V a r a h , J.M., On t h e S t a b i l i t y o f B o u n d a r y C o n d i t i o n s f o r Separable Difference Approximations to Parabolic E q u a t i o n s , SIAM J . Numer. A n a l . , 14 ( 1 9 7 7 ) 1114-1125. 19. V a r a h , J.M., S t a b i l i t y R e s t r i c t i o n s on S e c o n d O r d e r , T h r e e L e v e l F i n i t e D i f f e r e n c e Schemes f o r P a r a b o l i c S y s t e m s , SIAM J . Numer. A n a l . , 17 ( 1 9 8 0 ) 300-309.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The numerical solution of a system of diffusion/convection...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The numerical solution of a system of diffusion/convection equations describing coastal circulation Nicol, Thomas O. 1983
pdf
Page Metadata
Item Metadata
Title | The numerical solution of a system of diffusion/convection equations describing coastal circulation |
Creator |
Nicol, Thomas O. |
Date Issued | 1983 |
Description | A mathematical model describing nearshore ocean currents is examined. The motivation for the problem is discussed and a derivation of the model equations presented. The model equations consist of a pair of coupled, nonlinear, parabolic partial differential equations. A numerical method, the Crank-Nicolson finite difference scheme, is presented for solving the corresponding linear equations, and the convergence and stability of the method discussed. Methods for dealing with the nonlinear terms, and their effects on accuracy and stability, are examined. The initial and boundary conditions of the model equations present special problems, and we describe methods to solve them. A more efficient finite difference scheme for our equations, based on the Crank-Nicolson formula, is introduced, and its advantages discussed. A computer program incorporating these methods is developed to solve the model equations, and results for test data presented. We conclude with recommendations for additions to the program to model the currents more accurately. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0095808 |
URI | http://hdl.handle.net/2429/23990 |
Degree |
Master of Science - MSc |
Program |
Statistics |
Affiliation |
Science, Faculty of Statistics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1983_A6_7 N52.pdf [ 3.65MB ]
- Metadata
- JSON: 831-1.0095808.json
- JSON-LD: 831-1.0095808-ld.json
- RDF/XML (Pretty): 831-1.0095808-rdf.xml
- RDF/JSON: 831-1.0095808-rdf.json
- Turtle: 831-1.0095808-turtle.txt
- N-Triples: 831-1.0095808-rdf-ntriples.txt
- Original Record: 831-1.0095808-source.json
- Full Text
- 831-1.0095808-fulltext.txt
- Citation
- 831-1.0095808.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0095808/manifest