”'1 ‘ m" —— 7' 431-5- ~u.-rwu».~1‘-'I"-~~x. ':-‘" vv' - 1H . 401w #:V- . {(3. . ‘- , I ’ I . . 1' ' ‘P/.:uru""|\")”r’:1:r“’\h 2 .'. ' 31,191,. . Ir ' ' . ""3 "-‘ .. ,3” 1'11. 1111 Wu if?“ k K; ‘11-.” 3:33:11 1‘17. .131 1 “-1., fig: "" 37?. £11" . 'w , Add» 5"}:‘3'.’ “1'31 “‘b 314:1“ .‘fi‘i’t \f1- i -‘ - ' “7111?. ‘1‘ fig-1.. 1 ’. I ‘ _¥‘._Q:\\{V." 7:1}, ‘3’." hi‘ifi“:fi\‘-* *1 {_ “-‘r‘3;';r:.fr.“‘“” ' .. I . ‘ "~..1~..'I _ ‘ 3‘. ' .3 (if: 1 .. b ’- _"5“9~.§5 fl... , 1 ".A.._1.gh ‘ M ' ' ‘ S‘fiiafig‘fi ~ a. H- aggélqhfif . "ww. i 1 - '19?“ . . ' mm -»-_3- " “ ' 11.. fig?“ 11 1 fig 1 =»;:._§§% : 1w, 31-. - 1"”? 1 ~ - 1» a- 1 - -a «*r- .. 1:11 1 v t; "1‘ 353%» ) .g, ,1: r5 ‘1‘): ‘1 W}; 'x “If? $31"?! $14“. - ‘5' ‘ 'L ' l . ’ , ‘ ' .L‘Mépgg 97A; ‘0 [fit—'Wa-LV “K" '5 9}. 3:42.31‘ 1 $13.05? ‘11 1% $9,113.33 .thp :2! .' 1w .'. .‘.:h": :"(0- 3‘2 at: y" 1 wishe." 1,1 1}.“ "Qt-g. V .'. I h. ', ‘ n ‘ 1‘- ‘Q I _, “EMMAL'N ‘ £3531”- l LIBRARY Michigan State Unlvonlty llllllllllllllllll * v “€96 This is to certify that the thesis entitled NUMERICAL DERIVATION OF A MOMENTUM EQUATION OF HIGH FLOW RATE FLOWS IN POROUS MEDIA presented by Mikyoung Lee has been accepted towards fulfillment of the requirements for hastens—degreeinJlechgniLal Engineering at Major professor Date May 18, 1987 0-7839 MS U is an Afi‘irmau've Action/Equal Opponmu'ty [urination MSU ‘ LIBRARIES .—:—. RETURNING MATERIALS: Place in book drop to remove this checkout from your record. FINES will be charged if book is returned after the date stamped below. &P25fi% NUMERICAL DERIVATION OF A MOMENTUM EQUATION OF HIGH FLOW RATE FLOWS IN POROUS MEDIA By Mikyoung Lee AN ABSTRACT FOR A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Mechanical Engineering 1987 ABSTRACT NUMERICAL DERIVATION OF A MOMENTUM EQUATION OF HIGH FLOW RATE FLOWS IN POROUS MEDIA by Mikyoung Lee The objective of this thesis was to develop a relationship between the specific discharge or Darcian velocity and the pressure drop in a porous material when inertia effects are important. An array of circular cylinders with uniform dimension was assumed as a model of the porous medium. By considering the viscous incompressible flow in the narrow gap between circular cylinders, the dimensionless form of conservation equations were developed. These equations are solved in stream function form by quasi- linearizing, finite differencing, and applying successive over- relaxation to the resulting difference equations. The pressure drop across the cylinders is calculated from the stream function. It is found that the pressure drop is a quadratic function of specific discharge, in good agreement with experimental observations. ACKNOWLEDGEMENTS I wish to express my deepest appreciation to the members of my thesis committee, Professor Beck, Professor Bartholomew, and Professor Somerton, for their guidance and critical evaluation of this thesis. A special note of thanks goes to Professor Somerton for his patience, help and suggestion. I am deeply indebted to Professor Beck for sharing with him the wisdom needed to overcome several problems. His understanding and availability is appreciated. I also would like to express my appreciation to my friends who encouraged me a lot, especially Arafa Osman and Donhee Lee. Finally, I wish to express my greatest thanks to my father for his love, support, and understanding throughout the thesis work and all my years of study. TABLE OF CONTENTS page 1. LIST OF TABLES .................................. iii 2. LIST OF FIGURES ................................. iv 3 . NOMENCLATURE .................................... v CHAPTER 1 Introduction ....................................... 1 CHAPTER 2 Development of the Physical Model .................. 5 CHAPTER 3 Describing Equation and Boundary Condition ......... 20 CHAPTER 4 Method of Solution ................................. 29 CHAPTER 5 Results and Discussion ............................. 49 CHAPTER 6 Conclusions and Recommendations .................... 56 4. REFERENCES ...................................... 57 5. APPENDIX ........................................ 61 ii Tablel. Table2. Table3. Table4. List of Tables Title Representative values of porosity for various ...... substances Representative values of specific surface for ...... various substances Permeability models ................................ Difference equations ............................... iii Page 11 44 Figurel. Figure2. Figure3. Figure4. Figures. Figure6. Figure7. Figure8. Figure9. List of Figures Title Page The pressure drop-specific discharge data .......... 3 Example of a porous medium ........................ 6 Intergranular flow .................. ‘ .............. 13 Flow between Modelled parabola .................... l3 Flowchart ......................................... 47 Flowchart ......................................... 48 Effects of grid size on numericcal solution to .... 51 pressure drop for Reynolds number 0.5 Dimensionless pressure distribution ............... 53 Dimensionless pressure drop versus Reynolds ...... 55 number for empirical and numerical cases iv i,j J1’J'1HJ3 :1 t" 7‘: 7‘: q RR *- e, eA,Re . NOMENCLATURE surface area of the pores : diameter of grain the energy of flow per unit weight in x-direction gravity acceleration 1 piezometric head differential step size of n coordinate differential step size of 6 coordinate integer term indices : hydraulic gradient : hydraulic gradients along the axes : hydraulic conductivity : permeability 1 the average free distance between two grains : permeability 2 : pressure specific discharge Reynolds number total cross-sectional area the corresponding components of the velocity vector v X,Y,Z bulk volume volume of void space the volume of solids within vb the axes of a cartesian system of coordinates Greek letters atflifi19fi2 Subscripts X,y,Z numerical shape factors porosity density of the mixture dynamic viscosity kinematic viscosity specific surface specific area stream function transformed coordinates the relaxation parameter area the axes of a cartesian coordinate Superscript - : average value * : value with dimension, special index k : value of previous iteration vii CHAPTER 1 Introduction Many engineering processes involve the flow of a fluid through a porous medium. Examples include packed bed chemical reactors, degraded nuclear reactor cores, thermal and chemical pollution of aquifers, and secondary oil recovery in a petroleum reservoir. Historically, one may trace the orgins of an engineering concern for flows in porous media to the ancient Romans and their beautiful fountains as evidenced by the Trevi Fountain in Rome. These fountains‘ were once flow through devices which employed the pressurized water from the underground rivers of Rome. By standing the test of time these fountains demonstrate that the ancient Romans understood the engineering principles of flow through porous medium , but it was left to a French civil servant to express these principles in mathematical form. While designing the public water works for the village of Dijon, Darcy(l856) concluded that the specific discharge through a porous medium is directly proportional to the pressure drop across the medium. Darcy's law, as it has become known has, can be expressed mathematically for one dimensional flow in porous media as ' ____ ' q (1-1) where q is the specified discharge and K is the permeability, a property of the porous medium which will be discussed later. Though Darcy's law is a sound mathematical model for homogeneous flow through porous media there are a number of engineering applications when it is not appropriate. For non-homogeneous porous media Boussinesq(l904) and Irmay(l958) extended the linear equation of Darcy's. Similar extensions have been developed by Schneebeli(1957) and Ferrandon(l948) for non-isotropic media, by Muskat(l937) for compressible fluid flow, and by Anderson(1942) for solutions and adsorptive media. A detailed review of several modifications of Darcy's law is given by Muskat(l937). All of the physical processes listed above require some modification in Darcy's law. However the most serious breakdown of the mathematical model occurs when the flow rate is large. In Figure 1 the pressure drop-specific discharge data of Ahmed(l969) is plotted. It is seen that at small flow rates a linear relationship exists between pressure gradient and specific discharge. However, as the flow rate grows the data shows a marked deviation from the linear relationship of Darcy's law. In viewing this trend in his own experimental data, Forchheimer(l930) proposed a quadratic relationship between pressure gradient and specific discharge. The Forchheimer extension to Darcy's law can be expressed mathematically for one dimensional flow as '_"__q+_€l (1-2) where M has been called the Forchheimer coefficient and is dependent on the physical character of the porous medium. It is the theoretical basis of Eq(1-2) which is the main concern of this thesis. As it currently stands the Forchheimer extension may be considered semi- empirical and heuristic in nature. In this thesis first principles will be used to verify the form of Eq(1-2) for high flow rate flows in porous media, to identify the (DIMENSIONLESS) a , HYDRAULIC GRADIENT l . q, VELOCITY (cm/sec) FigL l I I I l l l l l l I l l l l I l I T I I I _ “ LEGEND 'i . . BROWNELL ' — I b FANCHER' '03 "_ C p — : o AHMED q }FORCHHElMER : '- o q KIRKHAM " h A v _. '02 _ o ALLEN ' }LINDQUIST I . BLAKE o MOBASHERI I ' 0 BROWNELL o SUNADA - l0 :— __ I —- _. IO" :— —_: '0'2 r -: Io-s _ .— D" J III] I 111' 1 i'ILI 14in L Jl '3 no-2 no" - no no2 physical process giving rise to the quadratic term, and to develop an expression for the Forchheimer coefficient in terms of the geometric parameters of porous media. This thesis continues with a discussion on the physical characteristics of porous media which will lead to the development of the physical model used in the thesis. The describing differential equations and boundary conditions of the physical model are then presented. The method of solution applied to the equations follows. The results of the solution are discussed and the thesis concludes by addressing the three points listed above. CHAPTER 2 Development of the Physical Model Before the physical model used in the thesis is presented a discussion of the physical characteristics of porous media is needed. The typical structure of a porous medium is shown in Fig.2. At any point in the interstitial space the flow obeys the Navier- Stokes equations. Due to the obvious complexity of the boundaries of the interstitial space a three-dimensional viscous flow in the tortuous channels of the medium arises. This occurs in spite of the apparent one-dimensional nature of the flow in terms of the specific discharge. It is clear from these observations there are two perspectives from which to view fluid flow in porous media: an interstitial perspective and a global or Darcian perspective. These two perspectives can be coupled through the following definition of the specific discharge as related to the interstitial velocity,’ q- _1_ J {”3 dA (2-1) where Ac is the total crossectional area, including pore space and solid, in the direction of the specific discharge and n is its unit normal vector. The preferred perspective for engineering analysis is the Darcian perspective. The two principle reasons involve the ease of measuring a specific discharge rather than an interstitial velocity and the difficulty in characterizing a porous medium in the interstitial Fig. 2 Example of a porous medium perspective. To fully characterize a porous medium from the interstitial perspective a complete description of the very complex boundaries of the interstitial space would be required. Clearly for a real porous medium this would be next to impossible. From the Darcian perspective average or bulk geometric properties are used to characterize the porous medium. It is generally agreed upon that three properties are necessary to fully characterize the medium. Perhaps the three most useful are the porosity,e, the specific surface, 2, and the tortuosity,r. The porosity is defined as the ratio of volume of void space to the bulk volume of the porous medium, 6 - (2-3) Some representative values of porosity for a variety of media are given in Table l. The specific surface is defined as the total interstitial surface area of the porous per unit bulk volume of the porous medium, 2 - (2-4) Some typical values of specific surface are given in Table 2. The tortuosity is defined as the square of the ratio of the distance traveled by a fluid element as it passes through porous sample to the length of the sample, TABLE 1 Representative Values of Porosity for various substances Substance (porosity in %) Berl saddles Raschig rings Wire crimps Black slate powder Silica powder Silica grains(grains only) Sand Granular crushed rock Soil Coal Concrete Leather Fibre glass Cigarette filters Hot-compacted copper powder Porosity range 68-83 56-65 68-76 57-66 37-49 65.4 37-50 44-45 43-54 2-12 2-7 56-59 88-93 17-49 9-34 Literature reference Carman,(1937) Ballard and Piret,(1950) Carman,(1937) Carman,(1937) Carman,(1937) Shapiro and Kolthoff(l948) Carman,(1937) Bernard and Wilhelm,(l950) Peerlkamp,(1948) Bond et al.,(l950) Berbeck,(l951) Mitton,(l945) Wiggins et al.,(1937) Corte,(1955) Arthur,(1956) TABLE 2 Representative values of Specific Surface for various substances Substance Specific surface range Literature reference (specific surface in cm- ) Berl saddles Raschig rings Wire Crimps Black slate powder Silica powder Catalyst Sand Leather Fibre glass 3.9-7.7 2.8-6.6 2.9*1o - 4.0*1o 3 3 7.0*1o - 8.9*10 3 3 6.8*10 - 8.9*10 S.6*105 2 2 1.s*1o - 2.2-10 a 4 l.2*10 - l.6*10 2 2 5.6*10 - 7.7*1O 1 Carman,(1937) Ballard and Piret,(l950) Carman,(1937) Carman,(1937) Carman,(1937) Spengler,(l936) Carman,(1937) Mitton et al.,(1945) Wiggins et al.,(1939) L 2 r - [ fe ] (2-5) L The two properties of momentum transport in porous medium, the permeability and the Forchheimer coefficient, must depend on these geometric properties. Formally, it may be written K - K(e,2,r) (2-6a) M - M(e,E,r) (2-6b) The exact form of the functions represented by Eqs(2-5) and (2-6) is difficult to specify for a real porous medium. A number of investigators have developed functional forms from their experimental data for the simple porous medium consisting of a bed of uniformed sized spheres. These forms are shown in Table3. It is important to note the variability in the models even with the simple geometry being considered. Several permeability models have been developed through theoretical derivations of Darcy's law. An excellent review of these models is given by Scheidegger (1960). Some of these derivations employ capillary-tube analogies, hydraulic radius theories by Terzaghi (1951) and a turbulent flow analogy by Yuhara (1954). These models have the following weak points: 1) they are mostly analogies, ii) their applicability have not been proved, or iii) they are based on the introduction of unspecified coefficients. 10 TABLE 3 Permeability Models Source K M E 1952 s d 2 3 d rgun( ) e p e p 2 150(1-6) l.75(l-e) 2 Schneebeli(1957) dp dp 1100 12 s 2 s 1 1 Carman(1953) e d e d ' P P 2 o 1 1 180(1-6) 2.87u ' (l-e) I (1958) 3 d 2 3 d rma y e p ‘ p 2 180(1-6) 0.5(1-6) 11 None of these models lead towards the Forchheimer extension. The physical source of the quadratic term in the Forchheimer extension has been argued considerably. One argument states that the term represents turbulence while another argument is that it represents inertia. As will be seen shortly the later argument seems to be correct. The only derivation of the Forchheimer term of any merit is that of Irmay(l958). Though not rigorous his derivation is provided here as background. Consider the packing of spherical particles shown in Fig.3. The Navier-Stokes equations for the flow of an incompressible Newtonian fluid through the interstitial space can be written as 0.: CI - - 2- + (u.V)u - - 1 VP + p V u (2-7) p n The continuity equation is V.u - 0 (2'8) Now introducing the following vector identity, (&.V)a - _1_ v<fi.a) - (6*(v*fi)) (2-9) 2 and substituting a a + _1_ v (a.a) - (&*(v*&)) - - _E_ 2- V P + v V u (2-10) 8 t 2 P 12 ////// Consider the overall flow in the x-direction and assume that the interstitial flow is steady and two-dimensional. By introducing the total energy of flow per unit weight of fluid as (2-11) its derivative in the X-direction may be written 6(gE) _ 1 a 2 2 2 (v ) - a (uv) + u a V + u [ a u + a “ ] (2-12) —— — — —-§ -—5 a x 2 ax ay 3 Y a x a Y This equation will be averaged vertical and denoting the vertical average with an overscore so that - 1 L F - Y F dy (2-13) L y 0 For the u velocity this will result in the specific discharge q - e G (2-14) where L J y u dy (2-15) 0 14 By reason of homogeneity and isotropy the Y-component of velocity must vanish upon averaging. That is, v V - 0 (2-16) Assuming that no correlation exists among the velocity components yields, - u - 0 (2-17) By applying continuity it follows that Bu 6v — - - — (2-188) 8X BY 2 2 6 u _ _ 6 v (2-18b) 6 X BXBY 2 2 8 u a v (2-l8c) 2 a X aan The energy will simply become the average energy on a vertical line given by E - (2-19) NI c1 00 13 E 15 where P is the average pressure along a vertical line. With these simplifications the momentum equation reduces to _ —-—-—a-- -r-- 6(gE) _ 1 a (v ) + V a u (2-20) —__2' a x 2 a x a Y The two remaining terms can be approximated from order of magnitude arguments which leads to the less than rigorous nature of the derivation. It is observed that in the interstitial space the u velocity will achieve a maximum as we travel along a vertical line so that we expect < 0 (2-21) It seems reasonable to have a “ (2-22) 2 1 or introducing the concept of permeability (2-23) 16 For the final term note that in the throat of the interstitial space v - 0 (2-24) Hence, in the converging zone 6 2 (v ) < 0 (2-25) 8X . in the diverging zone there will most probably be separation of flow at larger velocities, so that a 6X (v2) < 0 (2-26) is also valid there. Again an order of magnitude estimate can be employed using 2 _ 2 8(v ) - 11 (2-278) ax - a (2-27b) so that 2 - 2 31L L (2-28) 6 x a or introducing the Forchheimer term 17 ___ (2-29) 6 X With these expressions the momentum equation for porous media flow may be written as - - 2 _1_d_P_+_1_d(“)-- " q-i <2-30> p dX 2 dX K M But -2 - ‘1‘“) - 2 {i d_“_ (2-31) ax ax (2-32) dX The resulting equation is 2 +_”_q+Lq -0 (2-33) K M §2| %. which is the Forchheimer extension of Darcy's law. For the numerical derivation of a momentum equation for porous media flows presented in this thesis an appropriate model must be chosen. Though the spherical packing of Irmay's (1958) derivation 18 looks attractive, it is not strictly a two-dimensional model. Shayesteh(l984) used an array of cylinders as a porous media model to obtain a heat transfer coefficient relationship. He also showed that by neglecting inertia terms Darcy's law could be derived. This model,shown in Fig.4, will be used in the thesis. To be consistent with the description of a porous medium the porosity and specific surface of this model can be obtained as e - 1 - " (2-34) (1+L/d) z - " (2-35) d(1+L/d) These two expressions will allow one to represent the array of cylinders as a porous medium. 19 CHAPTER 3 Describing Equations and Boundary Conditions The describing equations for this problem are discussed considering the viscous incompressible flow in the narrow gap between two circular cylinders, see Fig.3. The active forces are those due to pressure and shear resulting from the fluid viscosity. The mathematical analysis is intended to define the wall equation and non-dimensional form of momentum equations. The continuity equation can be written - 0' (3-1) * * * 2* 2* * * p u 6 u + v 6 u _ _ 6 p + p 6 u + 6 u (3_2) * * 2 2 ax aY ax* ax* aY* and 'k * * 2* 2* * * p u 6 v + v 6 v _ _ 6 p + p 6 v + 6 v (3-3) 2 2 ax* aY* aY* ax* 311” The boundary conditions are : at the surface where Y* - H(X*) 20 u - O , v - 0 (3-4) ‘A’ 6 u * * - o , v - o (3-5) a Y aSX"'Q u* - Uo , v* - o (3-6) aSX"w u* - Uo , v* - o (3-7) * * To describe the boundary H(x ) consider a circle centered at Y - A+T. Beginning with the equation for a circle, 2 2 2 X + Y - A (3-8) * * and transforming X and Y to X and Y . X - X (3-9) 21 Y - Y*- T - A (3-10) Substituting and rearranging gives 2 2 2 Y* - 2Y*( A+T ) + x* - - T - 2TA (3-11) as the equation of the upper boundary. Perturbation and.Asylptotic Expansion.flethod To develop a simplified form of these equation the perturbation and asymptotic expansion method employed by Shayesteh(l984) is followed. Assume that there is some small parameter about which we can expand the dependent variables. It is found that o - T/A (3-12) will serve well as the small parameter. Thus the distance between the cylinders is considered small with respect to the diameter of the cylinder. Applying the following non- dimensionalization gives X - ----—- (3-13a) 22 Y Y - (3-13b) T 11: U. u - U0 (3'13C) * v v - v (3-13d) s . * p 3 13 p _ ( - e) P Order of magnitude arguments are used to determine a(Q), psand vs. Substituting the dimensionless variables the equation of the upper boundary becomes : 2 2 2 2 2 2 Y T - 2TY(A+T) + a A X - -T - 2TA (3-14) 2 dividing this equation by A , 2 2 2 2 2 2 Y O - 2Y(Q + O ) + a X - -¢ - 2O (3-15) The parameter a must be chosen to retain the physics of the problem. To begin let a - 1.0, then 23 2 2 2 2 2 X + Y O - 2Y(O + O ) - -O - 2O (3-16) Looking at terms of this equation in terms of the order of O 2 0(1) : X - 0 (3-17) ‘ 2 Since using X - O as the equation of the wall is meaningless,the 2 2 choice for a must be wrong. It is desired to have an a X in the order O terms so let a - O (3-18) Then 2 2 2 2 Y O - 2Y(O + O ) + OX - - O - 2O (3-19) The highest order terms in this equation are of O(O) which gives 2 -2Y + x - -2 (320) Then the equation for the wall becomes Y - H(X) - X2/2 + l - l/2(X2 + 2) (3-21) Note that this is equivalent to transforming the flow between the cylinders to that through the parabolic channel shown in Fig.4. 24 Using the functional form of a the non-dimensionalized continuity equation becomes U0 6 u vs 6 v , + . o (3'22) __172—— —-- - O A 6 X T 6 Y Rearranging, 6 u U0 1 6 v + . 0 (3—23) 1/2 _—_‘ ' 6 X v O 6 Y s 6 u 6 v The physics implies that and ___ should be of the same size to 6 X 6 Y retain the equality, hence vs 1 ___ z 0(1) z 1 (3-24) 1/2 Uo O 1/2 vs - O Uo (3-25) Now considering the dimensionless form of the momentum equation in * X - direction. 2 2 1/2 Uo 6 u Uo O 6 u 1 ps 6 p u + v - - . + ‘T72’ "“"" “‘ / O A 6 X T 6 Y p O A 6 X 25 u . + (3-26) The above equation can be written in the following form after being 2 multiplied by A O/VUo 1/2 1/2 2 AUOO 6 u 6 u pS A O 6 p 6 u 1 u +v - . + +_ u 6 X 6 Y p u Uo 6 X 6 X O 2 6 u (3-27) 2 6 Y ] 00 A where ReA - (3-28) V In words, Inertia force - Driving force + Viscous Force For a small Reynolds number the driving force must be balanced by the viscous force or 1/2 Ps A° z i (3-29) p u U0 O P U0 so that p - (3'30) 8 —'—3/7_ AO 26 Then in final form, A _— 2 2 1 2 Re O / u 6 u + v 6 u — - _l 6 p + a u + 1 6 u (3-31) O a X a x O a Y For the Y-momentum equation non-dimensionalizing gives 6 v ___ ___ __§.____ _ 6 X 1 2 Re O / u A 2 [ 6Y O6Y ax O6Y 2 2 + v 6 v _ _ l 6 p + 6 v + 1 6 v ] (3_32) "” ___§ In summary the dimemsionless governing equations for the flow between the cylinders are given as : continuity equation : ux + vY - 0 (3-33) momentum equation in X - direction : P X + “XX + uYY (3-34) O O 1/2 ReAO ( u ux + v uY) - - momentum equation in Y - direction : 5/2 2 ReAO ( u vx + v vY) - - pY + O vxx + O vYY (3-35) Assuming that the parameter O is small and taking only the highest order terms in the momentum equations gives 27 ReAO3/2( u ux + v uY) - - px + uYY (3-36) for the X-momentum equation. Since the analysis requires the retention of the inertia terms (otherwise one would expect to derive Darcy's law as Shayesteh(l984) did ) it is required that ReAO3/2 - 0(1) (3-37) or ReA - 0(2'3/2) (3-38) Then Re*( u ux + v uY) - -4px + uYY (3-39) and pY — 0 (3-40) will serve as the momentum equation where Re*- ReAOa/z - 0(1) (3-41) 28 CHAPTER 4 Method of Solution In the previous chapter the describing differential equations apprOpriate for the present problem were derived from the Navier- Stokes equations. To summarize the describing differential equations are ux + vY - o _ (4-1) Re* ('u uX + v uY ) - - px + uYY (4'2) pY - 0 (4-3) which must be solved subject to the boundary conditions ; at Y - 0 v - 0 , uY - 0 (4-4) at Y - h(X) u - O , v - 0 (4-5) 29 as X 4 - w u - l , v - 0 (4-6) as X 4 m u - 1 , v - 0 (4-7) Equations (4-1) through (4-3) represent a system of coupled, partial differential equations for the dependent variables u, v and p. To apply a numerical technique in solving these equations, manipulations are carried forth to yield a single partial differential equation of a single variable, the stream function. To begin, the pressure is eliminated by operating with g__ on Eqn (4-2) and noting pXY - (PY)x - 0 '(4-8) from the Y - momentum equation, then Re* ( uY ux + u uXY + vY uY + v uYY ) - uYYY (4-9) Now defining the stream function as u - ¢Y , v - - OX (4-10) 30 so that continuity is satisfied, substitution into Eq (4-9) gives * I”YYYY I Re ( Xx ¢YYY ' ¢Y l”YYX ) - 0 (4-11) which is a nonlinear partial differential equation for the stream function. The boundary conditions on the stream function are at Y - O AXO. at Y h(X) wXO. as X - w 2x0. as X m 0. 31 (4-12) (4-12) (4-13) (4-13) To handle the nonlinear terms, quasi-linearization is employed. The dependent variable is expanded by a Taylor's series about a previous various value ( from a previous iteration ). Thus, 0 px — ¢X + A Ox (4-l4a) 0 ¢YYY ' wYYY + A ¢YYY (“'IAb) 0 ¢Y - ¢Y + A wY (4-l4c) O ¢XYY ' ¢XYY I A I”XYY (“'14d) The superscript 0 denotes the value of the variable at the previous iteration. Substiting the expansions into Eq(4-10) and dropping higher order terms ( that is , products of A terms), * O O 0 0 * 0 O ¢YYYY+ Re ( ¢x MYYY + ¢YYY¢X ’ ¢Y wXYY ' I(’XYY””Y ) ' Re ( ¢YYY¢X ' o o ¢XYY¢Y 9 (4.15) The current X,Y coordinate system is somewhat burdensome, especially in light of finite differencing, due to the irregular boundary. To give a more regular computational domain the following coordinate transformation is employed. 32 1: - (4-16) h(X) 5 - x (4-17) 2 where, h(X) - x + 1 (4-18) 2 Therefore '1 - Y (4-19) T__ e /2 + 1 The derivatives can then be written 6 _ 6 n 6 + 6 5 6 (4-20a) 3—? a Y a 27 a Y a 5 - 1 a (4-20b) 2 ( 5 /2 + 1 ) 6 n a _ an a + as a (4-21.3) 5—): a x a 97 a x a e - 5 Y a + a (4-21b) 2 2 (6/2+1) an as 33 _ - 6 n 6 + 6 < 52/2 + 1) a n a s Define b - ( 52/2 + 1 )'1 %--enbwn+% a-bwfl 2 2 ¢XY ' ' 5 n b Nn" + b $5" ' 5 b V" 3 2 3 1bXYY - - 5 n b $000 + b $600 ' 2 6 b wan b2 UYY - Una b3 l(’YYY - lbdrm Substituting into the momentum equation gives 34 (4-21c) (4-22) (4-23a) (4-23b) (4-23c) (4-23d) (4-23e) (4-23f) (4-233) a * b3 o b3 o b3 o b3 2 Unnnn + Re ( we 1pmm ' wn Uénn + ¢nnn¢€ ' ¢€nn¢n + b‘wo w R * b8 o o b8 o o 2 b‘wo 2° 4 24 5 0 an) - e ( ¢nnn¢€ ' $0 $622 + 5 0 an ) ( ' ) The boundary conditions in the 6 and n coordinate system become at n - 0 $0" - O , $6 - 0 (4-25) at n - l O" - O , W5 - 0 (4-26) as 6 4 -w 1 w” - —B— . $5 - 6 0 (4-27) age-too w - 1 w - e n (4-28) n ‘g‘ ’ 6 A finite difference method will be used to solve Eq(4-24) subject to boundary conditions. Using a central difference approach gives 35 wi’J+l - Ibi’j_1 ¢ _ (4-29a) " 2 H 0 ¢ - ¢ $6 _ 1+1,j i-‘,j (4-29b) 2 H 5 4, -2)It +11) ¢ _ 1.3+1 1.1 1.3-1 (4-29c) an 2 H 2 ¢ l"1,.I+2 ‘ 2¢i.J+‘ + 2wi.J-‘ - ¢i»J'2 (4-29d) "WV 3 2 H n w _ ¢Y.J+2 ‘ “i1.3+1 * 6¢1.J ’ “ii-J-‘ + 21.3-2 (4-29e) nnnn 4 H n w _ (1+1,3+1 ' 2¢1+1,3 + l”1+j-1 '¢1-1,3+1 + 2X1-I,3‘ ¢i-1,j-1 £02 2 2H H n 6 (4-29f) Substitution and rearrangement gives a difference equation of the form 36 pi’j + C6 vi,j+1 + S Co $1-1.J-1 + C1 $1-1,J + C2 ¢1-1.j+1 + Ca P1’1-2 + C4 P1’3-1 + C5 i,j where, C? Ui,j+2 + Ca ¢i+1,j-1 + C9 ¢i+1,j+ C10 ¢i+1,j+1 ‘ (4-30) *0 R - e2" 2 H2 H n 6 'k *0 R R _ _ 3 $0 , 6 ¢ nnn H H 2 H 6 6 Re* w - 2 H2 H n 5 *0 b R - -6"? H‘ 2 H8 6 n *0 *0 * Re p Re P 2 b 5 Re p + 6 + £20 + 8 2 H 2 H H n n n 'k 6 b 4 e b Re ¢ - 4- 2 H H n n 37 (4-31a) (4-3lb) (4-3lc) (4-31d) (4-31e) (4-31f) Ce-“ ' + H H 2 H H n n n n 0 b Re* OE C7- ‘ + H 2 H 0 Re* O C8-- 0" 2 H H n E *0 Re* O Re O c, _ n + nun H H 2 H 6 f 0 Re* O c _ , n 10 2 2 H H n 6 S R * o o o o 2 b o o - - + 1.3 e ( $200 $6 $0 ¢€nn 5 $0 ¢nn) (4-313) (4-31h) (4-31i) (4313) (4-3lk) (4-311) The above equation will be valid for 3 S i s n-2 and 4 s j 5 m-3. To handle the other nodes the boundary conditions must be appl Numerically the boundary conditions as 5 4 i w are handled by ied. defining L585 half of the channel length. The correct size for L5 38 will have to be determined from numerical testing. Thus the boundary conditions as 6 4 i o become boundary conditions at 6 - i LE. Writing the boundary conditions in difference form gives at i - l 2 le-i-l ' M’l’jJ _ Lg + 1 (4-32a) 2 H 2 '7 1A ' '1’ 2.1 °.J _ - L n (4-32b) 2H 5 .1 6 at i - M 11» 1» L2 M'3+1 “'3'1 - 5 + 1 (4-33a) 2 H 2 ’7 ¢M+1.J - $M-1.J _ L5 "j (4 33b) 2 HE at j - 1 $1.0 - 2¢1,1 + ¢1,2 - o (4-34a) ¢i+1,1 + ¢i-‘,1 ' 0 (4’34b) 39 at j - N i”1,1~1+1 ' ¢i,N-1 ’ O (4'353) I"nan ’ $1-1,N ' 0 (“’35b) Since all of the boundary conditions are on the derivatives of the stream function we must choose to apply a boundary condition on the stream function itself. We choose O (6.0) - 0 (4-36) so that in difference form Next consider the i-l boundary where O1 3 — O1 1 + 2 H" ( LE/ 2 + 1 ) (4-38a) but O1 1 - O (4-38b) so that 40 O1 3 - 2 H"( L: / 2 + l ) (4-38c) Similarly 2 O1 5 - O1 3 + 2 H" ( L5/ 2 + 1 ) (4-38d) or substituting O, 5 - 4 H" ( LE/ 2 + 1 ) (4-38e) The trend is obvious so that 2 ¢l,N - ( N - 1 ) H" ( L5 / 2 + 1 ) (4-38f) but L" - ( N - 1 ) H" - 1 (4-38g) so that ¢,,N - L: / 2 + 1 (4-39) Equations (4-37) and (4-39) now replace Eqs. (4-34b) and (4-35b) at the boundaries. N1 Similarly it may be shown that at f - i L5 2 ib-(é /2+1)'I (4-40) Then Eqs (4-32a) and (4-33a) can be replaced with 2 O, 1 - "3 ( Lé / 2 + 1 ) (4-41) 2 OM . - "j ( LE / 2 + 1 ) (4-42) By way of Eqs(4-37), (4-39), (4-41) and (4-42) the stream function is known at all boundary nodes. For the inner nodes ( 2 s i 5 M-1 ) and 3 s j s N-2 ) Eq(4-30) will be appropriat. This leaves the nodes j-2 and j-N-l for 2 s i 5 M-1 left to be addressed. Applying Eqs.(4- 34a) and (4-35a) O1 0 - 2 A1,. - $1,. '(4-43a) (4-43b) I(LN-L1 ' I”LN-1 into Eq(4-30) gives the final two equations, 42 Co $1-1’1 + C1 $1-1’2 + C2 $1-1,3 + ( C4 + 2 C3 ) $1.1 + ( Cs ' Cs ) 21,2 + Cs Pi’s + C7 N1,4 + Ca P1+1’1 + C9 P1+1,2 + C1o ¢i+1,3 ' 51,2 (4-44) and Co Ui-1,N-2 + C1 Ui-1,N-1 + C2 ¢i-1,N + Cs P1,N_s + Ca wi,N-2 + ( C5 - C7 ) wi,N-1 + ( C6 + 2 C7 ) O1,N + C. ¢i+1,N-2 + C9 ¢i+1,N-1 + C10 S (4-45) ¢1+1,N ' i,N-1 which are valid for 2 s i 5 M-1. A summary of the appropriate difference equation at each node point is given in Table 4. Table 4 represents a system of algebraic equations for the discretized stream function, O(i,j). Solution to this system will require iteration since the source term, S1 1' contains derivatives of the stream function. Though this system could be solved in matrix form, this proves to be impractical due to the large size of the matrix and iterative nature of the problem. Each node equation could be written in such a form so that only the stream function at the node is at the current iteration while the remaining stream functions $1,j are evaluated at the previous iteration. For example, Eq.(4-30) may be rewritten, L13 O 0 O $1,J —E; ( Si,j ' C0 $1-1,J_1 ' C1 ¢1_1’j ' C2 ¢1-1,j+1 ' C3 0 O O O 0 O ¢i.j'2 ' C4 ¢i,j-1 ' C5 ¢l,j ' C6 wi,j+1 ' C1 ¢1’J+2 ' C8 ¢1+1,j-1 ' 0 0 0 Ca ¢1+1’j,1 ' C0 ¢1+193 ' C10 ¢i+1,j+1 ) where once again the superscript ° refers to the previous iteration. The method of successive over-relaxation (SOR) is used to iteratively solve this system of equations. Table 4 Difference Equations Node Points Equation 1-1 15 j s N (4-41) i-m 15 j S N (4-42) l< i 1 “ 1_____: D 1 Stop Lita [gData 2 J 4} Setup(ISTART)]l [ISolns I] so Je2,N2 Print DIP T 7 *** Warning! Solution Not Converged. *** (: Stop :) Figure 5 M7 FLDHCHARI 2 (Subroutine 1) (T Setup(ISTART) ) L. 1 - [ DATA 3 '1 s l«l,M PSI(I,1)-0.0 PSI(I,N) #4 VALUE 1 _— < Return :) FIDHCHARI h (Subroutine 3) 5 i313 / VALUE 3 <:_ Return f) Figure 6 “8 FIDUCHART 3 (Subroutine 2) VALUE 2 FIDUCHART 5 (Subroutine 4) (:iPrint if) s 1*1,M J+1,N <: Return VALUE 4 I l VALUE 5 3 CHAPTER 5 Results and Discussion Before presenting the results of this thesis, it is important to verify the correctness of the numerical solution. This is done by first considering the sum of the squares of the residuals. Since the difference equation which is solved is only an approximation of the differential equation an error in the solution will occur. The solution of the finite difference method is substituted into the differential equation and the residual, or difference that occurs, is evaluated at every node point. Each residual is then squared and the squares are summed. The sum of the squares of the residuals gives some measure of the correct choice of convergence criteria. For this thesis a point-wise convergence criteria is used by defining K K-l w - ¢ 11 11 (5-1) 13 x-1 and iterating until at every node point ‘ij is less than some prescribed value. It is found that by imposing cij < 10" (5-2) at every node gives a sum of the squares of the residuals between -17 -32 10 and 10 depending on Reynolds number. This gives some M9 confidence that the difference equation is being solved correctly and is a good approximation for the differential equation. The second aspect that needs to be considered with respect to the correctness of the numerical solution is the choice of grid spacing, A6 and An . Since the pressure drop will prove to bethe key parameter to use the numerical model in the investigation of porous media, its dependence on grid spacing is investigated. In Fig.7, the influence of 6 grid spacing on the pressure drop is shown. The pressure drop is determined at the same locations for each point and the ratio of An/A£ is also kept constant. It is seen that over a good range of A5 the pressure drop is relatively invariant with grid size. From this analysis grid spacings of A6 - 0.4 (5-3) An - 0.1 (5-4) were chosen and used for all of the remaining computer runs. As was stated above the key parameter for the porous media model is the pressure. The pressure is calculated from the momentum equations. In 5 and n coordinates the 6 momentum equation is written P b? b3 R * b2 b2 b3 3 (5 5 ' - - + - + - 6 60 n wan" e ( wflfl¢€ ¢€n¢n 6 ¢” ) ) while the n-momentum equation becomes P - o (5-6) 50 5.9... 0.0 m0 mmmEDZ moqozbm mOm dome mmDmmmma OH 20:340m #555232 20 MN_m 0E0 “_O Hounfim N**ANUV 0¢.0 00.0 0N0 070 00.0 F _ I r . L r 0 mHZmo NIH ._.< x I hop 1 Tom I. s n I I 1.0” _ _ b L . _ . 0.? c1080 EHDSSBHd 51 Then the pressure as a function of streamwise position (é-direction) is obtained by integrating Eq(S-S). A typical pressure profile is shown in Fig.8. The shape is what could be anticipated from the analytical results of Shayesteh(l984) when Re-O. Recalling that the primary thrust of this thesis is to determine the validity of the Forchheimer term for flows in porous media, it may be written - - 2 A P ” * P * (5-7) -* where u is the averaged or Darcian velocity. Using the non- dimensionalization of Chapter 3 gives - 2 ,2 A P - A é u + A 0 ReA u (5-8) 2L K m f Due to the non-dimensionalization u - 1.0 (5-9) so that Eq(5-8) could be written as A P - a + b Re (5'10) Hence a linear relationship between the numerically determined pressure drop and the Reynolds number would be expected if the Q 52 0N ZO_._.Dm_n/.:.m_0 mmDmmmmm mmmngmzméo 0.9m m_x momo mmsmmmmm mmmjzgmzmza 230mm >521. jwmwm210m Z<§mm_m 0.— 0.0 0.0 5.0 0.0 0.0 ¢.0 0.0 N0 L _ P L P _ _ L 50 0.0 c1080 BEDSSEHd SSEWNOISNHWIG 55 CHAPTER 6 Conclusions and Recommendations The following conclusions can be drawn from this thesis ; An array of cylinders proves to be an acceptable model as a porous medium The Forchheimer relation appears to be valid in the context of high velocity flows in porous media. The quadratic term of the Forchheimer relation is seen to represent the inertia forces of momentum transport. Several recommendations follow which have been generated from this thesis ; The slope and intercept of the pressure drop-Reynolds number curve need to be investigated as the physical configuration of the model is changed. In considering the slope, a differential equation for 6(AP)/6 Re in the limit of Re-0 should be developed to see the sensitivity the physical configuration has on the slope. The influence of turbulent flow in porous media can be studied using the present model by introducing a Reynolds stress term. 56 REFERENCES . Ahmed, N. and Sunada, D. K., ”Nonlinear flow in porous Media" J. of the Hydraulics Division, ASCE, vol.95, P.1847,1969 . Anderson, "Computational Fluid Mechanics and Heat transfer," New York, McGraw-Hall,l984 . Anderson,E.M.,"The Dynamics of faulting and dyke formation with applications to Britain," Oliver & Boyd Ltd., Edinburgh,l942 . Arthur,G., J. Inst. Metals 84, 327,1956 . Ballard, J. J., Piret, E. L., Ind. Eng. Chem. 42, 1088,1950 . Bear J.,"Dynamics of Fluids in Porous Media," New York, Elsevier, 1972 . Bernard, R. A., Wilhelm, R. H., Chem. Eng. Progr. 46, 223, 1950 . Boussinesq, J., J. Math. 10, 5 and 363, 1904 . Bond, R. L. et al., Fuel 29, No. 4, 83, 1950 57 10. ll. 12. l3. 14. 15. 16. 17. 18. Carman, P. C., "Fluid Flow through a Granular Bed," Trans. Inst. Chem. Eng. London 15, 150-156 1937 Corte, H.,"Z. Erzerg. Holzst.," lest., Papier und Pappe,9, 289, 1955 Darcy, H.,"Les fontaines publiques de la ville de Dijon," Dalmont, Paris,1856 Ergun, 8., Chem. Eng. Progr. 48, 89, 1952 Ferrandon, J., Genie Civil 125, 24, 1948 Forchheimer,”Formulas Transactions," American Geophysical Union 39(4) pp 702-707, 1930 Irmay 8.,"0n the theoretical derivation of Darcy", "Proc. Ass. Gen. Bruxelles," Ass. Int. Hydrol. (UGGI) 2, 179,1958 Merle C. Potter,” Mathematical Methods in the physical sciences," Prentice Hall, Inc., NJ, 1978 Muskat, M., ”The flow of homogeneous fluids through porous media," 1st en., 2nd prtg., publ. by Edwards, Ann Arbor, 1946 (763pp.), 1937 58 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. Mitton, R.G., J. Int. Soc. Leather Trades' Chem. 29, 255, 1945 Peerlkamp, P. K.,"Land bouwkund," Tijdschr. 60, 321, 1948 Piret, E.L. et al., Ind. Eng. Chem. 32, 861,1940 Scheidegger, Ph. D.,"The physics of flow through porous media," New York, The Macmillan Company,1960 Schneebeli, C., Pbl. Assoc. Int. Hydrologie (UGGI) 41, 10.: "Mem. Tr. Soc. Hydrotec.," France 2, 170, 1957 Shapiro, I., Kolthoff, M., J. Phys. Coll. Chem. 52, 1020, 1948 Shayesteh, M.,"A theoretical Model for the fluid-solid Heat transfer coefficient in porous media" M.S. Thesis, Louisiana State University,1984 Terzaghi, K., Proc. Am. Soc. Test. Mat. 45, 777.: "Theoretical soil Mechanics,” Chapman and Hall Ltd., London, 1951 Van Dyke,"Perturbation Methods in Fluid Mechanis," Annotated Edition, Stanford, Cal. 1975 Wiggins, E. J. et al., Canad. J. Res. B17, 318, 1939 59 29. Yuhara, K.,"Tikyubuturi Geophys. Inst.," Kyoto Univ. 9, No.2, 127,1954 60 APPENDIX '61 0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 0049 0050 0051 0052 0053 0054 0055 0056 0057 18-May-l987 10:4 lB-May-l987 10:3 PROGRAM FORCHHEIMER Ct*****t******t**.**i*ttfittttititttitttt*tfitti****i*tt*t*fl*ti*tit C‘ t C‘ DETERMINE PSI(Z,E) IN REGION I AM INTERESTED AND * C* THE RELATION BETWEEN PRESSURE DROP AND REYNOLDS NUMBER * Ct a Ct * Ctttt*tittttfltttttttitfittittittt*itttttttiittintfittittttitttttiit IMPLICIT REAL*8 (A-H.O-Z) REAL*8 LE,LZ COMMON/SOLS/PSI(100,50),PSIO(100,SO) COMMON/DERl/PSIOE(100,50),PSIOZ(100,50),PSIOZZ(100,SO) COMMON/DERZ/PSIOEE(100,50),PSIOZEE(100,50) COMMON/DER3/PSIOEEE(100,50),S(100,50),PREZ(100,50) COMMON/DER4/U(100,50),V(100,50) COMMON/DERS/PSIOZE(100,50),PREE(100,50) COMMON/DER6/PSIE(100,50),PSIZ(100,50) COMMON/GRID/M,N,M1,N1,N2,HE,HZ,LJ COMMON/PPARM/RE,BETHA,BIGDIF,DELP.LE,LZ COMMON/SQU/SQ(100,50) DIMENSION A3(100),H(100) Ct!***i***t*fltitttttfl*fi********t***********t*i*******t***tttfittti C* LEtLENGTH OF ETHA * C‘ LZ'LENGTH OF ZETHA * Cifittittfifitt**ttttfitfli******t**tiititflt**i*t*t****ti*ttt*titttflflt OPEN(UNIT=30,FILEI'MIKYOUNG.INP3',STATUSt'OLD') READ(30,21) NPTS 21 FORMAT(I2) DO 900 JJ'1,NPTS READ(30,22) RE.DELP,ISTRT,LSTJOB 22 FORMAT(F3.1,2X,F4.2,2X,12,2X,I) ctttttiiggtgttitttttttttttittt***ttt*******t****tt*tt**t*tttttt*tt C* IF ISTRT IS LESS THAN 0,CALCULATE INITIAL GUESS IN SETUP * C* IF ISTRT IS GREATER TRAN 0,USE PREVIOUS CALCULATION AT * C‘ LAST REYNOLD'S NUMBER INITIAL GUESSES. * C’ IF LSTJOB IS GREATER THAN 0, ONLY READING IN ONE SET OF * C‘ DATA. * * C***i**iflt****i*fiti*****.****t***********tfitfi**ttitttttttiitttiit OPEN(UNIT'35,FILEI'ANSWER3',STATUSI'NEW') BIGDIF-.01 M-22 N-lS LZ-3.0 LE-1.0 HZiZ.0*LZ/(M-l) HE-LE/(N-l) Nl-N-l N2-N-2 Ml-M-l 62 FORCHHEIMER 18-May-1987 10:4 lB-May-l987 13:3 0058 BETHA'.309 0059 WRITE(35,110) 0060 110 FORMAT(' ',3X,'DELP',5X,'BIGDIF',Gx,'REYNS 0') 0061 WRITE(35,120)DELP,BIGDIF,RE 0062 120 FORMAT(' ',1F7.2,3X,1F7.3,5X,1F7.3) 0063 WRITE(35,130) 0064 130 FORMAT(' ',3X,'HE',7X,'HZ',10X,'BETHA'I 0065 WRITE<35,140)HE,HZ,BETHA 0066 140 FORMAT(' ',1X,1F7.3,2X,1F7.3,5X,1F7.3) 0067 0068 0069 CALL SETUP(ISTRT) 0070 0071 - l8 CONTINUE 0072 0073 CALL SOLNS 0074 0075 DO 30 I81,H 0076 DO 30 J31,N 0077 0078 DIFIDABS((PSI(I,J)-PSIO(I,J))/PSIO(I,J)) 0079 IF(DIF.GT.1.D-4) GOTO 50 0080 30 CONTINUE 0081 0082 Ctttttt*ittittttttitttttitt*ttttt*tti*ttttttktttttttttttittt 0083 C* OUTPUT RESULT * 0084 Cttt*******ttttttttttttit*ttinttit*tttttttfittttttitt*ttttttt 0085 0086 CALL PRINT 0087 0088 GOTO 998 0089 50 CONTINUE 0090 0091 ctttit*ittittttttitn*ittt*******************tttt*tttittttttt 0092 C* NEW GUESSES BY S.O.R. * 0093 ctttttitttttttttt*tttttt*tttttitttttttttttitttttinttittttttt 0094 0095 DO 65 I'1,M 0096 DO 65 J31,N 0097 PSIO(I,J)-PSIO(I,J)-DELP*(PSIO(I,J)-PSI(I,J)) 0098 65 CONTINUE 0099 0100 CALL DERIV 0101 0102 LJ'LJ+1 0103 IF(LJ.GE.50000) GOTO 30 0104 0105 GOTO 18 0106 C1000 WRITE(35,90) 0107 C90 FORMAT(' ',"*‘ WARNING! SOLUTION NOT CONVERGED.**") 0108 0109 998 CONTINUE 0110 900 CONTINUE 0111 CLOSE(UNIT'30) 0112 CLOSE(UNIT'35) 0113 STOP 0114 END 63 lB-May-l987 10:4 18-May-l987 10:3 0001 SUBROUTINE SETUP(ISTRT) 883g Ciiifi***********t*i**********tt*ititfitfitit**i*********iti*tti***t 0004 C* SET STREAM FUNCTION AS INITIAL VALUES ‘ * 0005 Ciiitttttititti‘itt****t****t**ti******ititi***l**t**fi*********fi*i 0006 0007 IMPLICIT REAL*8 (A-H,O-Z) 0008 REAL'B LE,LZ 0009 0010 COMMON/SOLS/PSI(100.50),PSIO(100,50) 0011 COMMON/DERl/PSIOE(100,50),PSIOZ(100,50),PSIOZZ(100,50) 0012 COMMON/DERZ/PSIOEE(100,50),PSIOZEE(100,50) 0013 COMMON/DER3/PSIOEEE(100,50),S(100,50),PREZ(100,50) 0014 COMMON/DER4/U(100,50),V(100,50) 0015 COMMON/DERS/PSIOZE(100,50),PREE(100,50) 0016 COMMON/GRID/M,N,Ml,N1,N2,RE,HZ,LJ 0017 COMMON/PPARM/RE,BETHA,BIGDIF,DELP,LE,LZ 0018 COMMON/SQU/SQ(100,50) 0019 DIMENSION A3(100),H(100) 0020 0021 Al--2.0*((2.0)**(.5))/9.0/3.14 0022 Bl-Z.0*((2.0)**(.5))/3.0/3.l4 0023 0024 DO 5 I=1,M 0025 PSI(I,1)=0.01 0026 5 CONTINUE 0027 0028 DO 8 I=l,M 0029 PSI(I,N)=LE*(.5*((LZ)**2)+1.0) 0030 8 CONTINUE 0031 0032 IF (ISTRT.GT.0) GOTO 999 0033 0034 DO 6 I-1,M 0035 DO 6 J81,N 0036 PSIO(I.J)=(A1*(((J-l)*HE)**3)+Bl*(J-l)*HE)*45+0.01 0037 6 CONTINUE 0038 0039 DO 7I'l,M 0040 DO 7 J81,N 0041 PSIOE(I,J)'(-Bl*(((J-l)*HE)**2-1))*45 0042 7 CONTINUE 0043 0044 DO 11 I-1,M 0045 DO 11 J-1,N 0046 PSIOZ(I.J)=0.0 0047 11 CONTINUE 0048 0049 DO 9 I-1,M 0050 DO 9 J81,N 0051 PSIOEE(I,J)‘-Bl*(2.0*(J-l)*HE)*45 0052 9 CONTINUE 0053 0054 DO 10 I'I,M 0055 DO 10 J-1,N 0056 PSIOEEE(I,J)--2.0‘Bl*45 0057 10 CONTINUE 614 SETUP 0058 0059 0060 0061 0062 0063 0064 0065 0066 PROGRAM SECTIONS HOKDCDQO‘UIIFUNO FJH ENTRY POINTS 0-00000000 VARI 12 999 Name SCODE SLOCAL SOLS DERl DERZ DER3 DER4 DERS GRID PPARM SQU DO 12 I81,M DO 12 J81,N PSIOZEE(I,J)=0.0 CONTINUE CONTINUE RETURN END Total Space Allocated Address ABLES Address ** ARRAYS 2-00000000 R*8 Address Type Type Type R*8 10-00000018 R*8 AP-00000004@ 1*4 10-00000028 R*8 9-0000000C 1*4 Name SETUP Name A1 DELP ISTRT LZ N1 Name A3 2-00000320 R*8 H 8-00009C40 R*8 PREE Bytes 742 1608 80000 120000 80000 120000 80000 80000 40 48 40000 602438 Address Type 2-00000640 9-00000014 it 9-00000000 9-00000010 Bytes 800 800 40000 65 Attributes PIC PIC PIC PIC PIC PIC PIC PIC PIC PIC PIC CON CON OVR OVR OVR OVR OVR OVR OVR OVR OVR Name R*8 R*8 1*4 J 1*4 M 1*4 Dimensions (100) (100) (100, 50) REL REL REL REL REL REL REL REL REL REL REL LCL LCL GBL GBL GBL GBL GBL GBL GBL GBL SHR NOSHR SHR SHR SHR SHR SHR SHR SHR SHR SHR 10-00000008 9-0000001C 10-00000020 9-00000008 10-00000000 EXE NOEXE NOEXE NOEXE NOEXE NOEXE NOEXE NOEXE NOEXE NOEXE NOEXE Address 18-May-1987 10:4 18-May-l987 10:3 21503190211131?!)me T lB-May-19B7 10:4 18-May-1987 10:3 0001 SUBROUTINE DERIV 888g Cit*ittttfliitfitttfittt***t**t*******i**it****tfl*******fl****fi**ittt 0004 C* ALL DERIVATIVE OF STREAM FUNCTION ARE DERIVED BY * 0005 C* DIFFERENCE METHOD. * 0005 ctrtttittttttttttttttttttitttttttttttttittt*ttitttttttt*ttttttttt 0007 0008 IMPLICIT REAL*8(A-H,O-Z) 0009 . REAL*8 LE,LZ 0010 0011 COMMON/SOLS/PSI(100,50),PSIO(100.50) 0012 COMMON/DERl/PSIOE(100,50),PSIOZ(100,50),PSIOZZ(100,50) 0013 COMMON/DERZ/PSIOEE(100,50),PSIOZEE(100,50) 0014 COMMON/DER3/PSIOEEE(100,50),S(100,50),PREZ(100,50) 0015 COMMON/DER4/U(100,50),V(100,50) 0016 COMMON/DERS/PSIOZE(100.50).PREE(100,50) 0017 COMMON/GRID/M,N,M1,N1,N2,HE,HZ,LJ 0018 COMMON/PPARM/RE,BETHA,BIGDIF,DELP,LE,LZ 0019 DIMENSION A3(100),H(100) 0020 0021 DO 20 I-1,M 0022 0023 Z-(I-1)*HZ-LZ 0024 A3(I)-1.0/(.5*(2**2)+l.0) 0025 H(I)=(I-l)*HZ-LZ 0026 20 - CONTINUE 0027 0028 DO 21 I-2,M1 0029 DO 21 J'3,N2 0030 PSIOE(I,J)=(PSIO(I,(J+1))-PSIO(I,(J-1)))/2.0/HE 0031 21 CONTINUE 0032 0033 DO 22 I-2,M1 0034 DO 22 J-3,N2 0035 PSIOZ(I ,J)-(PSIO( (1+1) ,J)-PSIO((I-1),J) )/2.0/HZ 0036 22 CONTINUE 0037 0038 0039 DO 24 I-2,M1 0040 DO 24 J-3,N2 0041 PSIOZE(I,J)'(PSIO((I+1),(J+1))-PSIO((I-l),(J+1))-PSIO(( 0042 1 1+1),(J-1))+PSIO((I-1),(J-l)))/4.0/HE/HZ 0043 24 CONTINUE 0044 0045 DO 25 I-2,M1 0046 DO 25 J-3,N2 0047 PSIOEE(I,J)-(PSIO(I,(J+l))-2.0*PSIO(I,J)+PSIO(I,(J-1))1 0048 1 /((HE)**2) 0049 25 CONTINUE 0050 0051 DO 26 I-2,M1 0052 DO 26 J-3,N2 0053 PSIOEEE(I,J)'(PSIO(I,(J+2))-2.0*PSIO(I,(J+1))+2.0*PSIO 0054 1 (I,(J-1))-PSIO(I,(J-2)))/2.0/((HE)**3) 0055 26 CONTINUE 0056 0057 DO 27 I-2,M1 66 DERIV lB-May-1987 10:4 18-May-1987 10:3 0058 DO 27 J-3,N2 0059 PSIOZEE(I,J)=(PSIO((I+1),(J+1))-2.0*PSIO((I+1),J)+PSIO 0060 l ((1+1),(J-1))-PSIO((I-l).(J+1))+2.0*PSIO( 0061 2 (I-l),J)-PSIO((I~1),(J-1)))/2.0/((HE)**2)/HZ 0062 27 CONTINUE 0063 0064 DO 28 I-2,Ml 0065 PSIOE(I,1)=(-PSIO(I,3)+4*PSIO(I,2)-3*PSIO(I,1))/2/HE 0066 28 CONTINUE 0067 0068 DO 29 1'2,M1 0069 PSIOZ(I,1)'(PSIO((I+1),l)-PSIO((I-l),1))/2.0/HZ 0070 29 CONTINUE 0071 0072 0073 DO 31 I-2,M1 0074 PSIOZE(I,l)-(-(PSIO((I+1),1)-PSIO((I-l).l))/2/HZ+2*(PSIO((I+1). 0075 1 2)-PSIO((I-l),2))/HZ-3*(PSIC((I+l),3)-PSIO((I-l),3) 0076 2 )l2/HZ)/2/HE 0077 31 CONTINUE ‘ 0078 0079 DO 32 II2,M1 0080 PSIOEE(I,1)=(-PSIO(I,4)+4*PSIO(I,3)-5*PSIO(I,2)+2*PSIO(I,1)) 0081 1 /(HE)**2 0082 32 CONTINUE 0083 0084 DO 33 I-2,M1 0085 PSIOEEE(I,1)'(-3*PSIO(I,5)+l4*PSIO(I,4)-24*PSIO(I,3)+l8*PSIO 0086 1 (I,2)-5*PSIO(I,1))/2/(HE)‘*3 0087 33 CONTINUE 0088 0089 DO 34 I-2,M1 0090 PSIOZEE(I,1)-(-(PSIO((I+l),4)-PSIO((I-1),4))/2/HZ+2*(PSIO((I 0091 l +1),3)-PSIO((I-l),3))/HZ-2.5*(PSIO((I+l),2)-PSI 0092 2 O((I-l),2))/HZ+(PSIO((I+1),1)-PSIO((I-1),1))/HZ 0093 3 )/(HE)**2 0094 34 CONTINUE 0095 0096 DO 35 I-Z,M1 0097 PSIOE(I,N)'(3*PSIO(I,N)-4*PSIO(I,N1)+PSIO(I,N2))/2/HE 0098 35 CONTINUE 0099 0100 DO 36 I-2,M1 0101 PSIOZ(I,N)'(PSIO((I+1),N)-PSIO((I-l),N))/HZ/2.0 0102 36 CONTINUE 0103 0104 0105 DO 38 I-2,Ml 0106 PSIOZE(I,N)'(1.5*(PSIO((I+l).N)-PSIO((I-l),N))/HE-2*(PSIO((I-1 0107 1 ),N1)-PSIO((I-l),N1))/HZ+(PSIO((I+1),N2)-PSIO((I-l 0108 2 ).N2))/2/HZ)/2/HE 0109 38 CONTINUE 0110 0111 DO 39 I-Z,M1 0112 PSIOEE(I,N)-(2*PSIO(I,N)-5*PSIO(I,N1)+4*PSIO(I,N2)-PSIO(I,(N-3 0113 1 )))/(HE)**2 0114 39 CONTINUE 67 DERIV 0115 0116 0117 0118 0119 0120 0121 0122 0123 0124 0125 0126 0127 0128 0129 0130 0131 0132 0133 0134 0135 0136 0137 0138 0139 0140 0141 0142 0143 0144 0145 0146 0147 0148 0149 0150 0151 0152 0153 0154 0155 0156 0157 0158 0159 0160 0161 0162 0163 0164 0165 0166 0167 0168 0169 0170 0171 40 41 42 43 45 46 47 48 49. 50 52 UMP 18-May-1987 10:4 18-May-1987 10:3 DO 40 I32,M1 PSIOEEE(I,N)=(5*PSIO(I,N)-18*PSIO(I,Nl)+24*PSIO(I,N2)-14*PSIO (I,(N-3))+3*PSIO(I,(N—4)))/2/(HE)**3 CONTINUE DO 41 I32,M1 PSIOZEE(I,N)3 ((PSIO((I+1), N)- -PSIO((I- 1) ,N))/HZ- 2. 5*(PSIO((I+1 ) ,Nl)-PSIO((I- l) ,N1))/HZ+2*(PSIO((I+1), N2) -PSIO( (I- 1), N2))/HZ- (PSIO((I+1), (N- 3))- PSIO((I- 1L (N- 3 )))/2/HZ)/(HE)*'2 CONTINUE DO 42 1-2,M1 PSIOE(I,2)3(PSIO(I,3)—PSIO(I,1))/HE/2.0 CONTINUE DO 43 I32,M1 PSIOZ(I,2)3(PSIO((I+1),2)-PSIO((I-l),2))/2.0/HZ CONTINUE DO 45 132, Ml PSIOZE(I, 2)33(PSIO((I+1), 3)-4SIO((I-1L 3)-PSIO((I+1),1)+PSIO((I-l 1))/4. O/HE/HZ CONTINUE DO 46 I32,Ml PSIOEE(I,2)3(PSIO(I,3)-2.0*PSIO(I,2)+PSIO(1,1))/((HE)**2) CONTINUE DO 47 132, M1 PSIOEEE(I, 2)3 (- -3‘PSIO(I, 6)+14*PSIO(I, 5)- -24*PSIO(I, 4)+18*PSIO(I, 3 -5*PSIO(I, 2))/2/(HE)**3 CONTINUE W48 I32, M1 PSIOZEE(I, 2)3 3(PSIO((I+1), 3)-PSIO((I- 1L 3)- 2. 0*PSIO((I+1L 2)+ 2. 0*PSIO((I-1L 2)+PSIO((I+1L 1)- -PSIO((I- 1L 1)) /((HE)**2)/HZ/2. 0 CONTINUE DO 49 I32,M1 PSIOE(I,N1)3(PSIO(I,N)-PSIO(I,N2))/RE/2.0 CONTINUE DO 50 132, M1 PSIOZ(I, N1)3 3(PSIO((I+1) ,Nl)-PSIO((I- l) ,Nl))/HZ/2. 0 CONTINUE DO 52 132, M1 PSIOZE(I,N1)33(PSIO((1+1) ,N)-PSIO((I-l), N)-PSIO((I+1),N2)+PSIO((I N2))/4. O/HE/HZ CONTINUE DO 53 I32,M1 68 DERIV 0172 0173 0174 0175 0176 0177 0178 0179 0180 0181 0182 0183 0184 0185 0186 0187 0188 0189 0190 0191 0192 0193 0194 0195 0196 0197 0198 0199 0200 0201 0202 0203 0204 0205 0206 0207 0208 0209 0210 0211 0212 0213 0214 0215 0216 0217 0218 0219 0220 0221 0222 0223 0224 0225 0226 0227 0228 53 54 55 56 57 58 59 60 61 62 (AND-J 18-May-1987 10:4 18-May-1987 10:3 PSIOEE(I,N1)3(PSIO(I,N)-2.0*PSIO(I,N1)+PSIO(I,N2))/((HE)**2) CONTINUE DO 54 I32,M1 PSIOEEE(I,N1)3(5*PSIO(I,N1)-18*PSIO(I,N2)+24*PSIO(I,(N-3))-14*PS (I,(N-4))+3*PSIO(I,(N-5)))/2/(HE)**3 CONTINUE DO 55 I32,M1 PSIOZEE(I,N1)3(PSIO((I+1),N)-PSIO((I-l),N)-2.0*PSIO((I+1),N1) +2.0*PSIO((I-l),N1)+PSIO((I+1),N2)-PSIO((I-l),N2)) /2.0/HZ/((HE)**2) CONTINUE DO 56 J33,N2 PSIOE(1,J)=(PSIO(1,(J+1))-PSIO(1,(J-1)))/2.0/HE CONTINUE DO 57 J33,N2 PSIOZ(1,J)3(-PSIO(3,J)+4*PSIO(2,J)-3*PSIO(1,J))/2/HZ CONTINUE DO 1 J-3,N2 PSIOZE(1,J)3(-PSIO(3,(J+1))+4*PSIO(2,(J+1))-3*PSIO(1,(J+1))+ PSIO(3,(J-1))-4*PSIO(2.(J-1))+3*PSIO(1,(J-1))) /4/HE/HZ CONTINUE DO 58 J-3,N2 PSIOEE(1,J)3(PSIO(1,(J+1))-2.0*PSIO(1,J)+PSIO(1,(J-1)))/((HE)**2 CONTINUE DO 59 J33,N2 PSIOZEE(1,J)3((-PSIO(3,(J+1))+4*PSIO(2,(J+1))-3*PSIO(1,(J+1))) /2/HZ-(-PSIO(3,J)+4*PSIO(2,J)-3*PSIO(1,J))/HZ +(-PSIO(3,(J-1))+4*PSIO(2.(J-1))-3*PSIO(1,(J-1))) /2/HZ)/(HE)**2 CONTINUE DO 60 J33,N2 PSIOEEE(1,J)3(PSIO(1,(J*2))-2.0*PSIO(1,(J+1))+2.0*PSIO(1,(J-1)) -PSIO(1,(J-2)))/2.0/((HE)**3) CONTINUE DO 61 J33,N2 PSIOE(M,J)3(PSIO(M,(J31))-PSIO(M,(J-1)))/2.0/NE CONTINUE DO 62 J33,N2 PSIOZ(M,J)3(3*PSIO(M,J)-4*PSIO(M1,J)+PSIO((M-2),J))/2/HZ CONTINUE DO 2 J-3,N2 PSIOZE(M,J)3(3‘PSIO(M,(J+1))-4*PSIO(M1,(J+1))+PSIO((M-2),(J+l)) -3*PSIO(M,(J-1))+4*pSIo(M1,(J-1))-PSIO((M-2),(J-1)) )/4/HE/HZ CONTINUE D0 63 J-3.N2 69 DERIV 0229 0230 0231 0232 0233 0234 0235 0236 0237 0238 0239 0240 0241 0242 0243 0244 0245 0246 0247 0248 0249 0250 0251 0252 0253 0254 0255 0256 0257 0258 0259 0260 0261 0262 0263 0264 0265 0266 0267 0268 0269 0270 0271 0272 0273 0274 0275 0276 0277 0278 0279 0280 0281 0282 0283 0284 0285 63 64 65 wNH «gnaw H Nrd H «Jana H F‘ BJH 18-May-1987 10:4 18-May-l987 10:3 PSIOEE(M.J)3(PSIO(M,(J+1))-2.0*PSIO(M,J)+PSIO(M,(J-l)))/((HE)**2 CONTINUE DO 64 J33.N2 PSIOZEE(M,J)3((3*PSIO(M,(J+1))‘4*PSIO(M1,(J+1))+PSIO(("-2),(J+1) )/2/HZ-(3*PSIO(M,J)-4*PSIO(M1,J)+PSIO((M-Z),J))/HZ +(3'PSIO(M.(J-1))-4*PSIO(M1,(J‘1))*PSIO((M-2),(J-1) ))/2/HZ)/(HE)*‘2 CONTINUE DO 65 J33,N2 PSIOEEE(M,J)3(PSIO(M,(3+2))-2.0*PSIO(M,(J+1))+2.0*PSIO(M,(J-1)) -PSIO(M.(J-2)))/2.0/((HE)**3) CONTINUE PSIOE(1,1)3(-PSIO(1,3)+4*PSIO(1,2)-3*PSIO(1,1))/2/HE PSIOZ‘l.1)3(-PSIO(3.1)+4*PSIO(2.1)-3*PSIO(1.1))/2/HZ PSIOZE(1,1)3(-(-PSIO(3,3)+4*PSIO(3,2)-3*PSIO(3,1))+4‘(-PSIO(2.3) *4*PSIO(2,2)-3*PSIO(2,1))-3*(-PSIO(1,3)+4*PSIO(1,2) -3*PSIO(1,1)))/4/HZ/HE PSIOEE(1.1)3(-PSIO(1,4)+4*PSIO(1,3)-5*PSIO(1,2)+2*PSIO(1,1))/( HE *‘2 PSIOZEE(1,1)3('(-PSIO(3,4)+4*PSIO(2,4)-3*PSIO(1,4))/2/HZ-5*(- PSIO(3,2)+4*PSIO(2,2)-3*PSIO(1,2))/2/HZ+(-PSIO(3,1 )+4*PSIO(2,1)-3*PSIO(1,1))/HZ+2*(-PSIO(3,3)+4*PSIO( 2.3)-3*PSIO(1.3))/HZ)/(HE)**2 . PSIOEEE(1,1)3(-3*PSIO(1.5)+14*PSIO(1.4)-24*PSIO(1,3)*18*PSIO(1 ,2)-5*PSIO(1.1))/2/(HE)**3 PSIOE(1,2)3(PSIO(1,3)-PSIO(1,1))/NE/2.0 PSIOZ(1,2)3(-PSIO(3,2)+4*PSIO(2,2)-3*PSIO(1,2))/2/HZ PSIOZE(1,2)3(-(PSIO(3,3)-PSIO(3,1))+4*(PSIO(2,3)-PSIO(2,1))-3* (PSIO(1.3)-PSIO(1.1)))/4/HE/HZ PSIOEE(1.2)3(PSIO(1.3)-2*PSIO(1,2)+PSIO(1,1))/((HE)*‘2) PSIOEEE(1,2)3(-3*PSIO(1,6)+14'PSIO(1,5)-24*PSIO(1,4)+18*PSIO(1, 3)-5*PSIO(1.2))/2/(HE)**3 PSIOZEE(1,2)3((-PSIO(3,3)+4*PSIO(2.3)-3*PSIO(1.3))/2/HZ-(-PSIO (3,2)+4*PSIO(2.2)-3*PSIO(1.2))/HZ+(-PSIO(3.1)+4*PSI (2,1)-3*PSIO(1.1))/2/HZ)/(HE)‘*2 PSIOE(1,N)3(3*PSIO(1,N)‘4*PSIO(1,N1)+PSIO(1,N2))/2/HE PSIOZ(1,N)3(-PSIO(3,N)+4*PSIO(2,N)-3‘PSIO(1,N))/2/HZ PSIOZE(1,N)3(-(3*PSIO(3,N)-4*PSIO(3.N1)+PSIO(3,N2))/2/HE+4*(3* PSIO(2,N)-4‘PSIO(2,N1)+PSIO(2,N2))/2/HE-3'(3*PSIO( 1,N)-4*PSIO(1.N1)+PSIO(1,N2))/2/NE)/2/HZ PSIOEE(1,N)3(2’PSIO(1,N)-5‘PSIO(1,N1)+4*PSIO(1,N2)-PSIO(1.(N-3 )))/(HE)**2 PSIOEEE(1,N)3(5'PSIO(1,N)-18‘PSIO(1,N1)+24‘PSIO(1,N2)-14*PSIO(1 .(N-3))+3*PSIO(1.(N-4)))/2/(HE)'*3 PSIOZEE(1,N)3(2*(-PSIO(3,N)*4*PSIO(2,N)-3*PSIO(1,N))-5‘(-PSIO( 3,N1)+4*PSIO(2,N1)-3*PSIO(1,N1))+4'(-PSIO(3,N2)+4 ‘PSIO(2,N2)-3*PSIO(1,N2))-(-PSIO(3,(N-3))+4*PSIO( 2,(N-3))-3'PSIO(1.(N-3))))/2/HZ/(HE)*'2 PSIOE(1,N1)3(PSIO(1,N)-PSIO(1,N2))/HE/2.0 PSIOZ(1,N1)3(-PSIO(3.Nl)+4*PSIO(2,Nl)-3*PSIO(1,N1))/2/HZ PSIOZE(1,N1)3x-(PSIO(3.N)-PSIO(3,N2))/HE/2.0+2*(PSIO(2.N)-PSIO (2,N2))/HE-3'(PSIO(1,N)-PSIO(1,N2))/HE/2.0)/2/HZ 7O DERIV 18-May—1987 10:4 18-May-1987 10:3 0286 PSIOEB(1,N1)=(PSIO(l,N)-2*PSIO(1,N1)+PSIO(1,N2))/((HE)**2) 0287 PSIOEEE(1,N1)3(5*PSIO(1,N1)-18*PSIO(1,N2)+24*PSIO(1,(N-3))- 0288 1 14*PSIO(1,(N-4))+3*PSIO(1,(N-5)))/2/(HE)**3 0289 PSIOZEE(1,N1)= ((- -PSIO(3 N)+4*PSIO(2, N)- -3‘PSIO(1 N))/2/HZ- 0290 1 -(- -PSIO(3, N1)+4*PSIO(2 N1)- -3*PSIO(1 N1))/HZ+ 0291 2 (- PSIO(3, N2)+4*PSIO(2, N2)- -3*PSIO(1, N2))/2/HZ 0292 3 )/(HE)**2 0293 0294 PSIOE(M,1)=(-PSIO(M,3)+4*PSIO(H,2)-3*PSIO(M,1))/2/HE 0295 PSIOZ(M.1)=(3*PSIO(M.1)-4*PSIO(M1,1)+PSIO((M-2) l))/2/HZ 0296 PSIOZE(M 1)-(3*(-PSIO(M,3)+4*PSIO(M,2)-3*PSIO(M,1))-4*(-PS 0297 1 IO(M1 3)+4*PSIO(M1 2)- -3*PSIO(M1 1))+(- -PSIO((M- 0298 2 2), 3)+4*PSIO((M-2), 2)- -3*PSIO((M- 2) 1)))/4/HE/HZ 0299 PSIOEE(M,1)-(-PSIO(M 4)+4*PSIO(M, 3)- -5*PSIO(M, 2)+2*PSIO(M 1) 0300 1 )/(HE)**2 0301 PSIOEEE(M,1)=(-3*PSIO(M,S)+14*PSIO(M,4)-24*PSIO(M,3)+18*PSI 0302 2 O(M,2)-5*PSIO(M,1))/2/(HE)**3 0303 PSIOZEE(M,1)3(-(3*PSIO(M,4)-4*PSIO(M1 4)+PSIO((M- 2), 4))/2/HZ 0304 1 +2*(3*PSIO(M, 3)- -4*PSIO(M1, 3)+PSIO((M-2), 3))/HZ- 2 0305 2 .5*(3*PSIO(M, 2)- -4*PSIO(M1, 2)+PSIO((M-2), 2))/HZ+(3* 0306 3 PSIO(M,1)-4*PSIO(M1,1)+PSIO((M-2),1))/HZ)/(HE)**2 0307 0308 PSIOE(M,2)-(PSIO(M,3)-PSIO(M,1))/HE/2.0 0309 PSIOZ(M 2)- (3*PSIO(M, 2)- -4*PSIO(M1, 2)+PSIO((M- 2), 2))/2/HZ 0310 PSIOZE(M, 2)- (3*(PSIO(M 3)-PSIO(M, 1))- 4*(PSIO(M1, 3)- PSIO(M1, 1) 0311 1 )+PSIO((M-2), 3)- -PSIO((M-2) 1))/4/HE/HZ 0312 PSIOEE(M,2)-(PSIO(M, 3)-2*PSIO(M,2)+PSIO(M,1))/((HE)**2) 0313 PSIOEEE(M, 2)- (- -3*PSIO(M, 6)+14*PSIO(M, 5)- -24*PSIO(M, 4)+18*PSIO( 0314 1 M, 3)— —5*PSIO(M, 2))/2/(HE)**3 0315 PSIOZEE(M, 2)- -((3*PSIO(M 3)- -4*PSIO(M1. 3)+PSIO((M-2) 3))/2/HZ- 0316 1 (3*PSIO(M, 2)- 4*PSIO(M1, 2)+PSIO((M-2), 2))/HZ+( 0317 2 3*PSIO(M,1)-%*PSIO(M1,1)+PSIO((M-2),1))/2/HZ)/ 0318 3 (HE)**2 0319 0320 PSIOE(M, N1)3 (PSIo(M, N)-PSIO(M N2))/HE/2. 0 0321 PSIOZ(M, N1)3 (3*PSIO(M, Nl)-4*PSIO(M1, N1)+PSIO((M- 2) ,N1))/2/( 0322 1 H2) 0323 PSIOZE(M, N1)3 (3*(PSIO(M N)-PSIO(M, N2))- 4*(PSIO(M1, N)-PSIO(M 0324 1 1 N2))+PSIO((M-2), N)- 4SIO((M— 2) ,N2))/4/HE/HZ 0325 PSIOEE(M, N1)3 (PSIO(M, N)- -2*PSIO(M, N1)+PSIO(M, N2))/((HE)**2) 0326 PSIOEEE(M, N1)- (5*PSIO(M, N1)- -18*PSIO(M N2)+24*PSIO(M, (N- 3))- 0327 1 14*PSIO(M (N- -4))+3*PSIO(M, (N-5)))/2/(HE)**3 0328 PSIOZEE(M, N1)- ((3*PSIO(M, N)- 4*PSIO(M1, N)+PSIO((M-2) N))/2/HZ 0329 1 -(3*PSIO(M, N1)-4*PSIO(M1, N1)+PSIO((M-2), N1))/ 0330 2 Hz+(3*PSIO(M, N2)- 4*PSIO(M1, N2)+PSIO((M-2), N2) 0331 3 )/2/HZ)/(HE)**2 0332 0333 PSIOE(M, N)- (3*PSIO(M, N)- 4*PSIO(M, N1)+PSIO(M, N2))/2/HE 0334 PSIoz(M, N)- -(3*PSIO(M, N)- -4*PSIO(M1, N)+PSIO((M-2), N))/2/HZ 033s PSIOZE(M, N)-(3*(3*PSIO(M, N)- -4*PSIO(M.N1)+PSIO(M, N2))- 4*(3*PSIO 0336 1 (M1 N)- 4*PSIO(M1, N1)+PSIO(M1. N2))+3*PSIO((M-2), N)- 0337 2 4*PSIO((M-2) ,N1)+PSIO((M-2) ,N2))/4/HE/HZ 0338 PSIOEE(M, N)- (2*PSIO(M, N)- -5*PSIO(M, N1)+4*PSIO(M, N2)- PSIO(M, (N- 3) 0339 1 ))/(HE)**2 0340 PSIOEEE(M, N)- -(5*PSIO(M, N)- -18*PSIO(M N1)+24*PSIO(M, N2)-14*PSIO 0341 1 (M, (N-3))+3*PSIO(M, (N- 4)))/2/(HE)**3 0342 PSIOZEE(M, N)-((3*PSIO(M N)- 4*PSIO(M1 N)+PSIO((M-2) N))/HZ- 2. 5* 71 DERIV 0343 0344 0345 0346 0347 0348 0349 0350 0351 0352 0353 0354 0355 0356 0357 0358 0359 0360 77 DMNH 18—May-1987 10:4 18-May-1987 10:3 (3*PSIO(M,Nl)-4*PSIO(M1,Nl)+PSIO((M-Z),Nl))/HZ+ 2*(3*PSIO(M,N2)-4*PSIO(M1,(N-2))+PSIO((M-2),N2)) /HZ-(3*PSIO(M,(N-3))-4*PSIO(M1,(N-3))+PSIO((M-2), (N-3)))/HZ/2)/(HE)**2 DO 77 I=1,M DO 77 J-1,N ‘ PREZ(I.J)=((A3(I))**3)*PSIOEEB(I,J)-RE*(-H(I)*(A3(I))* *3*((PSIOE(I,J))**2)+(A3(I))**2*PSIOE(I,J)*PSIOZE (I,J)-(A3(I))‘*2‘PSIOZ(I,J)*PSIOEE(I,J)) CONTINUE RETURN END PROGRAM SECTIONS Name SCODE SLOCAL SOLS DERl DER2 DER3 DER4 DERS GRID PPARM ommqmmbwmo H Total Space Allocated ENTRY POINTS Address Type 0-00000000 VARIABLES Address Type 10-00000008 9-0000001C 9-00000024 9-00000004 ** Bytes Attributes 12408 PIC CON REL LCL SHR EXE R 1600 PIC CON REL LCL NOSHR NOEXE R 80000 PIC OVR REL GBL SHR NOEXE R 120000 PIC OVR REL GBL SHR NOEXE R 80000 PIC OVR REL GBL SHR NOEXE R 120000 PIC OVR REL GBL SHR NOEXE R 80000 PIC OVR REL GBL SHR NOEXE R 80000 PIC OVR REL GBL SHR NOEXE R 40 PIC OVR REL GBL SHR NOEXE R 48 PIC OVR REL GBL SHR NOEXE R 574096 Name DERIV Name Address Type Name Address T R*8 BETHA 10-00000010 R*8 BIGDIF 10-00000018 R*8 H2 ** 1*4 I ** 1*4 LJ 10-00000028 R*8 Lz 9-00000000 1*4 N 9-0000000C 1*4 N1 9-00000010 R*8 Z 72 0001 0002 0003 0004 0005 0006 0007 0008 0009 0010 0011 0012 0013 0014 0015 0016 0017 0018 0019 0020 0021 0022 0023 0024 0025 0026 0027 0028 0029 0030 0031 0032 0033 0034 0035 0036 0037 0038 0039 0040 0041 0042 0043 0044 0045 0046 0047 0048 0049 0050 0051 0052 0053 0054 0055 0056 0057 lB-May-1987 10:4 lB-May-1987 10:3 SUBROUTINE SOLNS C**t*************t*t**i***t****it*itttttfiitt***t******fi*t*******fl ct TO GET STREAM FUNCTION IN GIVEN CONDITION * Chittit*t**t*******i*******t***i**t****************i************* nnnn hJH Fl hJH (QRJH (AKJH IMPLICIT REAL*8 (A-H,o-Z) REAL*8 LB,LZ COMMON/SOLS/PSI(100,SO),PSIO(100,50) COMMON/DERl/PSIOE(100.50),PSIOZ(100,50),PSIOZZ(100,SO) COMMON/DERZ/PSIOEE(100,50),PSIOZEE(100,50) COMMON/DER3/PSIOEEE(100,50),S(100,SO),PREZ(100,50) COMMON/DER4/U(100,50),V(100,50) COMMON/DERS/PSIOZE(100,50),PREB(100,50) COMMON/GRID/M,N,M1,N1,N2,HB,HZ,LJ COMMON/PPARM/RE,BETHA,BIGDIF,DELP,LE,LZ COMMON/SQU/SQ(100,50) DIMENSION A3(100),H(100) DO 1 I'l,M DO 1 J81,N H(I)’(I-l)*HZ-LZ CONTINUE DO 10 J-2,N1 A3(l)-1.0/((LZ)**2/2.0+l.0) FF-(HE)**4*HZ*(.00001) co-(RB*PSIOE(1,J)/2.0/(HE)**2/Hz)*FF c1=-2.o*co-(RE*PSIOEEE(1,J)/2.0/HZ)*FF czaco c3-(A3(1)/(HE)**4-RE*Psxoz(1.J)/2.0/(HE)**3)*FF C4-(-4.0*A3(l)/(HE)**4+RE*PSIOZ(l,J)/(HB)**3+RE* PSIOZEE(1,J)/2.0/HE+RE*PSIOE(1,J)*2.0*H(1)*A3 (1)/(HE)**2)*PF C5-(6.0*A3(1)/(HE)**4-4*RE*H(1)*A3(1)*PSIOE(1,J) /(HE)**2)*FF C6-(-4.0*A3(1)/(ns)**4-RB*PSIoz(1,J)/(HE)**3-R£* PSIOZEB(1,J)/2.0/HE+2*RE*H(1)*A3(1)*PSIOE(1,J) /(HB)**2)*PF C7-(A3(1)/(HE)**4+RE*PSIOZ(l.J)/2.0/(HE)**3)*FF C8--C0 c9--c1 Clo-C8 S(1,J)-RE*(PSIOEEE(1,J)*PSIOZ(1,J)-PSIOE(1,J)*PSIOZEE(1,J) +2*H(1)*A3(1)*PSIOE(1,J)*PSIOEE(1,J))*FF PSI(1,J)'(S(1,J)-2.0*(C0*(J-2)+C1*(J-l)+C2*J)*HE*HZ*LZ '(C7-C3-C4)*HE*2.0/A3(1)-(C4+C6)*PSIO(1.(J+l) )-(C8+C0)*PSIO(2,(J-l))-(C9+Cl)'PSIO(2,J)-(C10+C2)*PSIO (2.(J+l)))/(C5+C3+C7) SQ(l,J)-(PSI(1,J)*(C5+C3+C7)+2*(C0*(J-2)+C1*(J-1)+C2*J)*HE*HZ*LZ +(C7-C3-C4)*HE*2/A3(1)+(C4+C6)*PSIO(1,(J+1))+(C8+C0)* PSIO(2,(J-1))+(C9+C1)*PSIO(2,J)+(C10+C2)*PSIO(2,(J+l)) -S(1.J))**(2) 73 SOLNS 18-May-1987 10:4 lS-May-1987 10:3 0058 10 CONTINUE 0059 0060 A3(M)-A3(l) 0061 DO 11 J82,Nl 0062 CO-(RE*PSIOE(M,J)/2.0/(HE)**2/HZ)*FF 0063 Cl'-2.0*C0-(RE*PSIOEEE(M,J)/2.0/HZ)*FF 0064 C2'C0 0065 C3-(A3(l)/(HE)*‘4-RE*PSIOZ(M,J)/2.0/(HE)**3)*FF 0066 C4=(-4.0*A3(l)/(HB)**4+RB*PSIOZ(M,J)/(HE)**3+RE*PSIOZEB(M,J) 0067 l /2.0/HE+2*RE*H(M)*A3(H)*PSIOE(M,J)/(HE)**Z)*FF 0068 C5t(6.0*A3(l)/(HE)**4-4*RE*H(M)*A3(M)*PSIOE(M,J)/(HE)**2)*FF 0069 C6‘(-4.0*A3(l)/(HE)**4-RE*PSIOZ(M,J)/(HE)**3-RE*PSIOZEE(M,J) 0070 l /2.0/HE+2*RE*R(M)*A3(M)*PSIOE(M,J)/(HE)**2)*FF 0071 C7-(A3(1)/(HE)**4+RE*PSIOZ(M,J)/2.0/(HE)**3)*FF 0072 C8--C0 0073 C9--Cl 0074 Clo-C8 007$ S(M,J)'RE*(PSIOEEE(M,J)‘PSIOZ(M,J)-PSIOE(M,J)*PSIOZEE(M,J) 0076 1 +2*H(M)*A3(M)*PSIOE(M,J)‘PSIOEE(M,J))‘FF 0077 0078 PSI(M,J)'(S(M,J)-(C8*(J-2)+C9*(J-1)+C10*J)*2.0*HE*HZ*LZ 0079 l -(C7-C3-C4)*HE*2.0/A3(l)-(C0+C8)*PSIO((M-l),(J- 0080 2 l))-(C1+C9)*PSIO((M-l),J)-(C2+C10)*PSIO((M-l),(J+l))- 0081 3 (C4+C6)*PSIO(M,(J+1)))/(C5+C3+C7) 0082 C SQ(M,J)I(PSI(M,J)*(C5+C3+C7)+(C8*(J-2)+C9*(J-1)+C10*J)*2‘HE*HZ* 0083 C 1 LZ+(C7-C3-C4)*HE*2/A3(1)+(C0+C8)‘PSIO((M-l),(J-l))+(C1 0084 C 2 +C9)*PSIO((M-l),J)+(C2+C10)*PSIO((M-l),(J+l))+(C4+C6)* 0085 C 3 PSIO(M,(J+1))-S(M,J))‘*(2) 0086 11 CONTINUE 0087 0088 0089 DO 15 I-2,M1 0090 Z'(I-1)*Hz-LZ 0091 A3(I)-1.0/(.5*(z**2)+1.0) 0092 0093 C0-(RE*PSIOE(I,2)/2.0/(HE)**2/HZ)*FF 0094 Cl--2.0*C0-(RE*PSIOEEE(I,2)/2.0/HZ)*FF 0095 CZ‘CO 0096 C3-(A3(I)/(HE)**4-RE*PSIOZ(I,2)/2.0/(HE)**3)*FF 0097 C4-(-4.0*A3(I)/(HE)**4+RE*PSIOZ(I,2)/(HE)**3+RE*PSIOZBE(I,2) 0098 1 /2.0/HE+2*RE*H(I)*A3(I)*PSIOE(I,2)/(HE)**2)*FF 0099 C5-(6.0*A3(I)/(HE)**4-4*RE*H(I)*A3(I)*PSIOE(I,2)/(HE)**2)*FF 0100 C6'(-4.0*A3(I)/(HE)**4-RE*PSIOZ(I,2)/(HE)**3-RE*PSIOZEE(1,2) 0101 l /2.0/HB+2*RE*H(I)*A3(I)‘PSIOE(I,2)/(HE)**2)*FF 0102 C7-(A3(I)/(HE)**4+RE*PSIOZ(I,2)/2.0/(HE)**3)*FF 0103 - C8--C0 0104 C9--Cl 0105 Clo-CB 0106 S(I,2)'RE*(PSIOEEE(I,2)*PSIOZ(I,2)-PSIOE(I,2)*PSIOZEE(I,2) 0107 1 +2*H(I)'A3(I)*PSIOE(I,2)*PSIOEE(I,2))*FF 0108 0109 PSI(I,2)'(S(I,2)-C1*PSIO((I-l),2)-C2*PSIO((I-1),3) 0110 1 -(C4+2.0*C3)*PSIO(I,1)-C6*PSIO(I,3)-C7*PSIO(I,4) 0111 2 -(C8+C0)*PSIO((I+1),1)-C9*PSIO((I+1),2)-C10*PSIO((I 0112 3 +1).3))/(C5-C3) 0113 C SQ(I,2)'(PSI(I,2)*(C5-C3)+C1*PSIO((I-l),2)+C2*PSIO((I-1),3)+ 0114 C 1 (C4+2*C3)*PSIO(I,1)*C6*PSIO(I,3)+C7*PSIO(I,4)+(C8+C0 7L: SOLNS 0115 0116 0117 0118 0119 0120 0121 0122 0123 0124 0125 0126 0127 0128 0129 0130 0131 0132 0133 0134 0135 0136 0137 0138 0139 0140 0141 0142 0143 0144 0145 0146 0147 0148 0149 0150 0151 0152 0153 0154 0155 0156 0157 0158 0159 0160 0161 0162 0163 0164 - 0165 0166 0167 0168 0169 0170 0171 Hnnnn GOO 18-May-1987 10:4 18-May-1987 10:3 )*PSIO((I+1),1)*C9*PSIO((I+l),2)+C10'PSIO((I+1).3)- S(I,J))**(2) CONTINUE um DO 16 I'Z,M1 C0=(RE*PSIOE(I,Nl)/2.0/(HE)**2/HZ)*FF C18-2.0*C0-(RE*PSIOEEE(I,N1)/2.0/HZ)*FF C2-C0 C3-(A3(I)/(HE)**4-RE*PSIOZ(I,Nl)/(HE)**3/2.0)*FF C4'(-4.0*A3(I)/(HE)**4+RE*PSIOZ(I,Nl)/(KE)**3+RE*PSIOZEE(I,Nl) 1 /2.0/HE+2*RE*H(I)*A3(I)*PSIOE(I,N1)/(HE)**2)*FF C5'(6.0*A3(I)/(HE)**4-4*RE*H(I)*A3(I)*PSIOE(I,N1)/(HE)**2)*FF C6=(-4.0*A3(I)/(HE)‘*4-RE*PSIOZ(I,N1)/(HE)**3-RE*PSIOZEE(I,Nl) 1 /2.0/HE+2*RE*H(I)*A3(I)*PSIOE(I,N1)/(HE)**2)*FF C7=(A3(I)/(HE)**4+RE*PSIOZ(I,N1)/2.0/(HE)**3)*FF C83-C0 C9I-C1 C10'C8 S(I,Nl)-RE*(PSIOEEE(I,N1)*PSIOZ(I,N1)-PSIOE(I,Nl)*PSIOZEE(I,N1) l +2*H(I)‘A3(I)*PSIOE(I.N1)'PSIOEE(I,N1))‘FF PSI(I,N1)-(S(I,Nl)-C0*PSIO((I-l),N2)-C1*PSIO((I-1).N1) -C3*PSIO(I,(N-3))-C4*PSIO(I,(N2))-C6*PSIO(I,N) -C8*PSIO((1+1),(N2))-C9*PSIO((I+1),N1)-(C10+C2)*PSIO( (I+1),N))/(CS+C7) SQ(I,N1)'(PSI(I,N1)*(C5+C7)+C0*PSIO((I-l),N2)+C1*PSIO((I-1),N1) 1 +C3*PSIO(I,(N-3))+C4*PSIO(I,N2)+C6*PSIO(I,N)+C8*PSIO( 2 (1+1),N2)+C9*PSIO((I+1),N1)+(C10+C2)*PSIO((I+1),N)- 3 S(I,Nl))**(2) UNI-4 CONTINUE DO 20 I-2,M1 DO 20 J'3,N2 2-(1-1)*HZ-LZ A3(I)-l.0/(.5*(z**2)+1.0) C0'(RE*PSIOE(I,J)/2.0/(HE)**2/HZ)*FF Clt-2.0*C0-(RE*PSIOEEE(I,J)/2.0/HZ)*FF C2-C0 C3-(A3(I)/(HE)**4-RE*PSIOZ(I,J)/(53)**3/2.0)*PF C4'(-4.0*A3(I)/(HE)**4+RE*PSIOZ(I,J)/(KE)**3+RE*PSIOZEE(I,J) 1 /2.0/HE+2*RE*H(I)*A3(I)*PSIOE(I,J)/(HE)**2)*FF C5-(6.0*A3(I)/(HE)**4-4*RE*H(I)*A3(I)*PSIOE(I,J)/(HE)**2)*FF C6-(-4.0*A3(I)/(HE)**4-RE*PSIOZ(I,J)/(HE)**3-RE*PSIOZEE(I,J) 1 /2.0/HE+2*RE*H(I)*A3(I)*PSIOE(I,J)/(HE)**2)*FF C7-(A3(I)/(HE)**4+RE*PSIOZ(I.J)/2.0/(HE)**3)*FF C8--C0 C9--C1 Clo-C8 S(I,J)-RE*(PSIOEEE(I,J)*PSIOZ(I,J)-PSIOE(I,J)*PSIOZEE(I,J) +2*H(I)*A3(I)*PSIOE(I,J)*PSIOEE(I,J))*FF PSI(I,J)-(S(I,J)-C0*PSIO((I-1).(J-l))-C1*PSIO((I-l),J)-C2* PSIO((I-1),(J+l))-C3*PSIO(I,(J-2))-C4*PSIO(I,(J-1)) -C6*PSIO(I,(J+1))-C7*PSIO(I,(J+2))-C8*PSIO((1+1),(J -1))-C9*PSIO((I+1),J)‘C10*PSIO((I+1),(J+1)))/C5 SQ(I,J)-(PSI(I,J)*C5+C0*PSIO((I-l).(J-1))+C1*PSIO((I-1),J)+C 1 2*PSIO((I-1),(J+1))+C3*PSIO(I,(J-2))+C4*PSIO(I,(J-1) 2 )+C6*PSIO(I,(J+1))+C7*PSIO(I,(J+2))+C8*PSIO((1+1),(J “NH H 75 SOLNS 0172 0173 0174 0175 0176 0177 0178 0179 0180 0181 0182 O OQOOODNOO O H -1))+C9*PSIO((I+1),J)+C10*PSIO((I+1),(J+1))-S(I,J))* *(2) CONTINUE SUM80.0 DO 101 I'l,M DO 101 J'1,N SUM'SUM+SQ(I,J) CONTINUE WRITE(20,*) RE,LJ,SUM RETURN END PROGRAM SECTIONS Name SCODE SLOCAL SOLS DERl DER2 DER3 DER4 DERS GRID PPARM SQU HOOQQO‘U‘FNMO Hid Total Space Allocated ENTRY POINTS Address Type 0-00000000 VARIABLES Address 10-00000008 2-00000698 2-00000670 2-00000690 9-0000001C 9-00000024 9-00000004 ** Bytes 5392 1704 80000 120000 80000 120000 80000 80000 40 48 40000 607184 Attributes PIC CON REL PIC CON REL PIC OVR REL PIC OVR REL PIC OVR REL PIC OVR REL PIC OVR REL PIC OVR REL PIC OVR REL PIC OVR REL PIC OVR REL Name SOLNS Type Name Address Type Name R*8 BETHA 10-00000010 R*8 BIGDIF R*8 C10 2-00000658 R*8 C2 R*8 C5 2-00000678 R*8 C6 R*8 C9 10-00000018 R*8 DEL? 3*8 Hz 2-000006A0 1*4 I 1*4 LJ 10-00000028 R*8 Lz 1*4 N 9-0000000C 1*4 N1 R*8 z 76 LCL LCL GBL GBL GBL GBL GBL GBL GBL GBL SHR EXE NOSHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE SHR NOEXE Address 2-00000648 2-00000660 2-00000680 2-00000640 2-000006A4 9-00000000 9-00000010 lB-May-1987 10:4 lB-May-1987 10:3 wwmwwmmmmmm T lB-May-1987 10:4 lS-May-1987 10:4 0001 SUBROUTINE PRINT 8883 Ctiittittttfltiiittttittittitittttttttttt*tittt*ttttttttttttittttt 0004 C‘ PRINT OUT THE RESULTS WE CONCERNED. * 0005 C*********fi**************tfi****t*ttti‘k‘kttt*tttfiifiiiifiiiifitii‘tiiit 0006 0007 IMPLICIT REAL*8 (A-H,O-Z) 0008 REAL*8 LZ,LE 0009 COMMON/SOLS/PSI(100,50),PSIO(100,50) 0010 COMMON/DERl/PSIOE(100,50),PSIOZ(100,50),PSIOZZ(100,50) 0011 COMMON/DERZ/PSIOEE(100,50),PSIOZEE(100,50) 0012 COMMON/DER3/PSIOEEE(100,50),S(100,50),PREZ(100,50) 0013 COMMON/DER4/U(100,50),V(100,50) . 0014 COMMON/DERS/PSIOZE(100,50),PREE(100,50) 0015 COMMON/DERS/PSIE(100,50),PSIZ(100,50) 0016 COMMON/GRID/M,N,M1,N1,N2,HE,HZ,LJ 0017 COMMON/PPARM/RE,BETHA,BIGDIF,DELP,LE,LZ 0018 DIMENSION A3(100),H(100),PRE(100,50),PRE1(100,50) 0019 0020 C WRITE(35,200) 0021 C200 FORMAT(' ',1X,'RE',2X,'I',2X,'J',2X,'PSI(I,J)') 0022 0023 C DO 20 I-2,M1 0024 C DO 20 J'3,N2 0025 0026 C PSIE(I,J)-(PSI(I,(J+1))-PSI(I,(J-1)))/2.0/HE 0027 C20 CONTINUE 0028 0029 C DO 21 1'2,M1 0030 C DO 21 J-3,N2 0031 C PSIZ(I,J)'(PSI((I+1),J)-PSI((I-1),J))/2.0/HZ 0032 C21 CONTINUE 0033 0034 C DO 22 I-2,M1 0035 C PSIE(I,1)-(2*PSI(I,2)-1.5*PSI(I,1)-.5*PSI(1,3))/HE 0036 C22 CONTINUE 0037 0038 C DO 23 I-2,M1 0039 C PSIZ(I,1)'(PSI((I+1),1)-PSI((I-l),1))/2.0/HZ 0040 C23 CONTINUE 0041 0042 C DO 24 I-2,M1 0043 C PSIE(I,N)'(1.5*PSI(I,N)-2*PSI(I,N1)+.5*PSI(I,N2))/HE 0044 C24 CONTINUE 0045 0046 C DO 25 I-2,M1 0047 C PSIZ(I,N)'(PSI((I+1),N)-PSI((I-1),N))/HZ/2.0 0048 C25 CONTINUE 0049 0050 C DO 26 I-2,M1 0051 C PSIE(I,2)-(PSI(I,3)-PSI(1,1)HHS/2.0 0052 C26 CONTINUE 0053 0054 C DO 27 I-2,M1 0055 C PSIZ(I,2)'(PSI((I+1),2)-PSI((I-1),2))/2.0/HZ 0056 C27 CONTINUE 0057 77 PRINT lB-May-1987 10:4 ‘ lB-May-l987 10:4 0058 C DO 28 I=2,M1 0059 C PSIE(I,N1)=(PSI(I,N)-PSI(I,N2))/HE/2.0 0060 C28 CONTINUE 0061 0062 C DO 29 I-2,M1 0063 C PSIZ(I,N1)'(PSI((I+1),N1)-PSI((I-1),N1))/HZ/2.0 0064 C29 CONTINUE 0065 0066 C ~ DO 30 J83,N2 0067 C PSIE(1,J)'(PSI(1,(J+l))-PSI(1,(J-l)))/2.0/HE 0068 C30 CONTINUE 0069 0070 C DO 31 J83,N2 0071 C PSIZ(1,J)'(PSI(2,J)-PSI(1,J))/HZ 0072 C31 CONTINUE 0073 0074 C DO 32 J-3,N2 0075 C PSIE(M,J)=(PSI(M,(J+1))-PSI(M,(J-1)))/2.0/HE 0076 C32 CONTINUE 0077 0078 C DO 33 J83,N2 0079 C PSIZ(M,J)=(PSI(M,J)-PSI((M-l),J))/HZ 0080 C33 CONTINUE 0081 0082 C PSIE(1,1)'(PSI(1,2)-PSI(1,1))/HE 0083 C - PSIZ(1,1)=(PSI(2,1)-PSI(1,1))/HZ 0084 C PSIE(1,2)'(PSI(l,3)-PSI(1,1))/HE/2.0 0085 C PSIZ(1,2)‘(PSI(2,2)-PSI(1,2))/HZ 0086 C PSIE(1,N)'(PSI(l,N)-PSI(1,N1))/HE 0087 C PSIZ(1,N)-(PSI(2,N)-PSI(1,N))/KZ 0088 C PSIE(1,N1)=(PSI(1,N)-PSI(1,N2))/HE/2.0 0089 C PSIZ(1,N1)'(PSI(2,N1)~PSI(1,N1))/HZ 0090 C PSIE(M,1)'(PSI(M,2)-PSI(M,1))/HE 0091 C PSIZ(M,1)-(PSI(M,1)-PSI(M1,1))/HZ 0092 C PSIE(M,2)'(PSI(M,3)-PSI(M,1))/HE/2.0 0093 C PSIZ(M,2)'(PSI(M,2)-PSI(M1,2))/HZ 0094 C PSIE(M,N1)-(PSI(M,N)-PSI(M,N2))/HE/2.0 0095 C PSIZ(M.N1)'(PSI(M,N1)-PSI(M1,N1))/HZ 0096 C PSIE(M,N)'(PSI(M,N)-PSI(M,N1))/HE 0097 C PSIZ(M,N)'(PSI(M,N)-PSI(M1,N))/HZ 0098 0099 0100 0101 DO 6 I=I,M 0102 Z-(I-1)*HZ-LZ 0103 A3(I)-l/(Z**2/2.0+1.0) 0104 6 CONTINUE 0105 0106 C DO 7 I-1,M 0107 C DO 7 J-1,N 0108 C U(I,J)-A3(I)*PSIE(I,J) 0109 C7 CONTINUE 0110 0111 C DO 8 I-1,M 0112 C DO 8 J-1,N 0113 C V(I,J)-(I-1)*HZ*(J-1)*HE*A3(I)*PSIE(I,J)-PSIZ(I,J) 0114 C8 CONTINUE 78 PRINT lB-May-1987 10:4 lB-May-1987 10:4 0115 0116 C DO 9 I=1,M 0117 C DO 9 J81,N 0118 C WRITE(35,210) RE,I,J,PSI(I,J) 0119 C210 FORMAT(' ',1P3.1,1X,I2,2X,I2,3X,Pl7.6) 0120 C WRITE(35,211) I,J,PREZ(I,J) 0121 C211 FORMAT(' ',12,2X,12,2X,1F7.3) 0122 C9 CONTINUE 0123 0124 C DO 99 I-1,M 0125 C DO 99 J-1,N1 0126 C PRE(I,J)-HE*(PREE(I,J)+PREE(I,(J+1)))/2.0 0127 C99 CONTINUE 0128 0129 DO 100 I=1,M1 0130 DO 100 J'1,N 0131 PRE1(I,N)=HZ*(PREZ(I,N)+PREZ((I+1),N))/2.0 0132 100 CONTINUE 0133 0134 C DO 101 J'1.N1 0135 C SUM1-SUM1+PRE(1,J) 0136 C101 CONTINUE 0137 C DO 102 J-1,N1 0138 C SUMZ-SUM2+PRE(M,J) 0139 C102 CONTINUE 0140 0141 C DO 79 I-1,H 0142 C SUM1-0.0 0143 C D0 101 J-l,N1 0144 C SUMI'SUM1+PRE(I,J) 0145 C101 CONTINUE 0146 OPEN(UNIT-31,FILE-'PRESSURE',STATUS-‘NEW') 0147 C WRITE(31,13) 0148 C13 FORMAT(' ',‘RE',5X,'I',3X,'PRESSURE') 0149 C WRITE(31,14) RE,I,SUM1 0150 C14 FORMAT(' ',1F3.1,3X,12,2X,1F17.9) 0151 C79 CONTINUE 0152 DO 103 I-l,M1 0153 DO 103 J'1,N 0154 WRITE(35,104) I,PRE1(I,N) 0155 104 FORMAT(' '.12,2X,1F17.6) 0156 103 CONTINUE 0157 0158 C DO 105 I'l,M1 0159 C WRITE(31,106)RE,PSIOE(I,1),PSIOZEE(I,1),PSIOZE(I,1),PSIOEE(I,1) 0160 C106 FORMAT(' '.1F3.1,2X,1F17.6,2X,1F7.3,2X,1F7.3,2X,1F7.3) 0161 C105 CONTINUE 0162 RETURN 0163 END 79 “11111011111