TRANSIENT TWO-PHASE FLOW THROUGH POROUS MEDIA Thasis for flm 5.3?er of Ph. D. Mtéfiififd‘é Sfifffi Pifiik’fifismf Kuang-mfing 5:351 WM: THESIS This is to certify that the thesis entitled Transient Two-Phase Flow Through Porous Media presented by Kuang-ming Lin has been accepted towards fulfillment of the requirements for fig)— degree in _Me_chanic s 46% 5% ajor professdr Date November 25, 196k 0-169 ROOa’A {ESE OHY ROGM ABSTRACT TRANSIENT TWO-PHASE FLOW THROUGH POROUS MEDIA by Kuang-ming Lin This thesis presents an analytical investigation by which the movement of the interface for transient two- phase flow through porous media can be determined. Two approaches, one-dimensional flow and two-dimensional flow, are considered in the analysis of the problem. In the one-dimensional flow, the governing equations are solved simply by a finite-difference method. It is found that the results are not very satisfactory near the outflow seepage face. Therefore, two-dimensional flow is emphasized. The essential idea utilized in two-dimensional flow is to treat the interface as a distribution of sources and to apply the concept of Green's function to the governing differential equation in obtaining the solution. For purposes of illustration, a case of con- fined finite aquifer of rectangular shape is considered. The two interacting fluids are assumed incompressible and immiscible. A Control Data 3600 computer is used for all the numerical computation. TRANSIENT TWO-PHASE FLOW THROUGH POROUS MEDIA BY Kuang—ming Lin A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Metallurgy, Mechanics, and Materials Science 1964 ACKNOWLEDGEMENTS The author wishes to express his gratitude to Dr. George E. Mase for his constant encouragement and advice throughout the author's graduate work at Michigan State University. Also, thanks are due Dr. Harold R. Henry for suggesting the subject of this dissertation and for his continuous guidance during its preparation. Dr. Lawrence E. Malvern, Dr. Joachim E. Lay, and Dr. J. Sutherland Frame are also to be thanked for their serving on the Guidance Committee and for reviewing this thesis. For the computational work of this thesis, the author is indebted to the Computer Center and its staff for their co-operation and generous allowance of the computer time which made the execution of the tremendous amount of computation in this thesis possible. Finally, the author would like to thank his wife, Pimyau, for her patience and encouraging words during the entire period of his graduate study. ii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . ii LIST OF TABLES O O O O O O O 0 0 O O O O O O O O 0 iv LIST OF FIGURES O O O 0 O O 0 O 0 0 O O O O O O O 0 V LIST OF APPENDICES O O O O O O O O 0 O o O O 0 O 0 Vi Chapter I 0 INTRODUCTION 0 O O 0 O 0 O O 0 0 O O O O O l A. Statement of the Problem . . . . . l B. General Theory—-Darcy's Law. . . . 4 C. Review of Literature . . . . . . . 12 II o ONE-DIMENSIONAL ANALYSIS 0 o o o o o o o o 15 A. Derivation of the Equations. . . . 15 B. Numerical Example. . . . . . . . . 20 C. Finite-difference Solution Procedure. . . . . . . . . . . . . 23 D. Results and Discussion . . . . . . 24 III 0 TWO-DINIENSIONAL ANALYSIS 9 o o o o o o o o 34 A. Formulation of the Boundary-value PrOblem. O 0 O 0 0 O o O O O O O O 34 B. Solution Procedure with Green's Function . . . . . . . . . . . . . 40 C. Numerical Example. . . . . . . . . 45 D. Results and Discussion . . . . . . 46 IV. CONCLUSION O O O O O O O O O O O O O O O O 6 2 BIBLIOGRAPHY O O O O O O O O O 0 O O O O O O O O O 64 APPENDICES 0 O O O O O O 0 0 O O O O O O O 0 O O O 66 iii Table LIST OF TABLES Tabulation for Figures 5 and 6 Tabulation for Figure 7. . . . Tabulation for Figure 8. . . . Values of PSI along Initial Interface for (Tabulation for y/d = 0.1 to 0.8 Figure 13) . . . . . . . . Movement of Interface Position—~Case 4 (Tabulation for Figure 14) Values of PSI for Initial Interface Position y = 0.810 to 0.982 (Tabulation of Figure 15). Movement of Interface Position-—Case (Tabulation of Figure 16). iv 0 0 O Page 31 32 33 52 54 55 61 Figure l. 3. 4. 6. 7. 8. 9. 10. ll. 12. 13. 14. 15. 16. Geometry LIST OF FIGURES of the Problem. . . . Generalized Darcy's Experiment . Graph of Flow Geometry for One-dimensional Equation (1). . . . . Analysis . . . . . . . . . Plot of h/d vs. x/d. . . . . . Movement Movement Movement Geometry Boundary Boundary Boundary Plot of Movement of Interface—-Case l. of Interface-~Case 2. of Interface-—Case 3. O of Flow and Boundary Conditions Conditions in Terms of yl. . . Conditions in Terms of 3p Conditions for G. . . I st.y/d...... of Interface--Case 4. W, VS. Y/do o o o o o o o o 0 Movement of Interface--Case S. O Page 16 27 28 29 30 35 39 4O 41 48 49 SO 51 LIST OF APPENDICES Appendix page I 0 LIST OF SYMBOLS O O 6 O O o 0 O 0 O O O O O 66 II. FORMULAS USED IN NUMERICAL APPROXIMATION . 68 A. Newton's Divided-difference Interpolation Formula. . . . . . . 68 B. Simpson's One-third Rule for Numerical Integration. . . . . . . 69 C. Numerical Differentiation FormUlaS o O 0 O O 0 O O O 6 O O O 69 III. FLOW CHART FOR COMPUTATION . . . . . . . . 70 IV. SAMPLE FORTRAN PROGRAMS. . . . . . . . . . 72 PROGRAM PSI 0 O O O 0 9 o 0 O O 0 O O 0 73 PROGRAM NEW POSN . . . . . . . . . . . 78 V. GREEN'S FUNCTION FOR LAPLACE'S EQUATION AND POISSON'S EQUATION . . . . . . . . 79 A. Definition . . . . . . . . . . . . 79 B. Application to the Problem of Chapter III 0 O O 0 0 O O O 0 O O 0 83 vi CHAPTER I INTRODUCTION A. Statement of the Problem The object of this thesis is to investigate analytically the unsteady flow patterns which exist in porous media when two adjacent fluids of different densities are in a non-equilibrium state. The primary consideration is given to the determination of the movement of the interface separating the fluids. The analytical treatment of such a problem is generally very difficult, and only a handful of solutions of special cases are available. Some of the more important of these with respect to this investigation are described in part C of this chapter and are also listed in the Bibliography. Practical engineering situations in which the results of this investigation may be applied occur in the problems of sea-water encroachment into fresh-water aquifers and in the area of petroleum recovery techniques. In this thesis, the physical arrangement in which the problem is formulated consists of a rectangular porous medium, confined in the horizontal direction, together with two adjacent liquid reservoirs as indicated schema- tically by Figure 1. For the problems to be considered 1 2 the horizontal dimension of the rectangle is actually some twenty times the vertical dimension, so that the inter- faces and the flow are much more nearly horizontal than they appear in Figure l. ___.L___4 reservoir 1 ljll/////////[//1//_1Ill/llljllllllllL/IIJ/ -"r .L ' ‘ ' I I H " I, ‘ d — ~ ‘ \ ’ 4 U, ‘I ~ - I v .I .‘ f ' x .s r 7 ' - ‘- W's - ‘- . _' ‘ .: ‘. " ‘ ’ . flnal _ ~ - ~‘ ~ ~ ‘ , ' interface ‘ ' . . _ - ‘ — _- n'i‘tiaflt. - '~ -‘- " , ' nterface~ l, ‘ ’ , , ~”—Zoné 22-, , _ - o - - ‘ ‘- 4 ' I "' J .4 - / I I I ‘ ' ~- ‘5 reservoir 2 l7/’//?//// // //7,I/x/I/UIA////III/III/ III //////////////// Figure l.--Geometry of the Problem When the piezometric head of one of the liquid reservoirs is suddenly changed by a certain amount and then kept at a constant value, the flow of liquid through the porous medium is unsteady and the interface moves toward a new equilibrium position. Such a condition is observed fre— quently in nature when a heavy rainfall occurs or during a drought when the heavy withdrawal of ground water supply takes place. It is the purpose of this thesis to apply an analytical treatment to the problem of the establishment of flow through a porous medium and especially the motion 3 of the interface under such conditions during the transient state. Darcy's law is a fundamental law upon which the mathematical theory of flow through porous media is formulated. The investigation as reported herein analyzes the problem from two different assumptions, one-dimensional flow and two-dimensional flow. The former is presented in Chapter II and the latter in Chapter III. In the one— dimensional analysis, the directions of the velocities are assumed to be in the x—direction only. Since the hori- zontal dimension is twenty times the vertical in Figure 1, this would seem to be a reasonable assumption over most of the region. The resulting governing differential equations are replaced by the corresponding finitewdifference equations, and these are solved numerically. The results thus obtained show that this approach is not completely satisfactory, and the need for tw0ndimensional analysis is apparent. The inaccuracy of the one-dimensional approach stems from the fact that the velocities near the out-flow seepage face change their directions very rapidly, thus making the one-dimensional flow assumption invalid. The emphasis, therefore, will be placed on the two- dimensional analysis, which describes more realistically the physical situation. The one-dimensional analysis is not totally without merit, however, since it does apply over the major part of the region and also serves as a guide to the two-dimensional analysis. 4 The essential idea employed in the two~dimensional analysis is to consider the interface as a distribution of sources and to apply the concept of a Green's function to the governing differential equations in obtaining a solution. The resulting solution is in the form of a double infinite series. With the aid of a modern high speed computer, such as the Control Data 3600 system at Michigan State University, sufficient convergence can be obtained in a reasonable amount of time to yield a solution within the range required by engineering accuracy. B. General Theory-—Darcy's Law Flow through a porous medium, like any other type of flow, obeys Newton's second law of motion, which states that "forces must be exerted on a fluid to change either the direction or magnitude of the fluid velocity." When a fluid flows through a porous medium, the velocity of a fluid element changes rapidly from point to point along its tortuous flow path. The forces which produce these changes in velocity vary rapidly from point to point. However, in a naturally porous material the porous structure and hence the multitude of flow paths have a random character. It is reasonable, therefore, to suppose that the random variations in flow patterns for any partin cular fluid element are uniformly distributed. Also the variations in magnitude of velocity can be expected to be distributed uniformly with mean zero. Thus, for steady laminar flow the lateral forces associated with the 5 microscopic random variations in velocity can be expected to average to zero over any macroscopic volume. However, the inertial forces in the direction of flow will not average to zero and hence will only be negligible for low flow rates. Fortunately, the flow in most cases of prac- tical interest is of the slow laminar type and Darcy's law, which is presented below, applies. The mathematical theory of flow through porous media is always formulated with Darcy's law being taken as the fundamental law of flow. In the middle of the 19th century, Henry Darcy, a French engineer, discovered through experiments a law governing the flow of water through filter beds. This law expressed in vectorial form is ‘5’ -=grad (0' (l) or -V - (2) is the velocity potential, and h is the piezometric head. For the meaning of‘$ and y, refer to Figure 2. Here, the liquid of specific weight 'r is flowing with a flow rate of Q (dimension L3/T) through a tube which is filled with ‘A complete list of symbols with their definitions is given in Appendix I. a porous medium of length L. 6 It is seen that .2 r and y represent the pressure head and the elevation, respectively. 1? E \ _ Ah=h2-h : 7 N":7 1 El 3 : p2 3‘ - : 7 ' _ Q .1 KW -’..\ h .. 2 y2 From dimensional analysis, Figure 2.-~Generalized Darcy's Experiment conductivity k can be expressed as k = cd 21 A o the hydraulic (3) The dimensionless factor c combines the effects of porosity, range and distribution of grain sizes, and shape of grains as well as their orientation and packing, while d, the mean grain size, is representative of the average pore size. The hydraulic conductivity k depends on the properties c and d of the medium as well as on the specific weight 7 and viscosity Au of the fluid. 7 The product K = cd2 is typical of the medium alone and is called the "permeability." Darcy's law is invariant with respect to the direction of flow in the earth's gravitational field. This can be proved by the flow in the apparatus sketched in Figure 2. This property was not recognized immediately, since Darcy performed his experiment in a vertical pipe. Figure 3 illustrates the vector character of Equation (1) with the velocity potential of Equation (2) for the case where )’= f>g is independent of position, as in a homogeneous incompressible fluid. Figure 3.--Graph of Equation (1) Since g‘7y = -g3, where g is the acceleration of l‘ gravity and j is a unit vector in the +y-direction, it follows from Equation (2) that for a homogeneous incompressible fluid " 1% grad ‘9 = 93 3%,- grad p- (4) 8 Equation (4) states that the force exerted upon the unit mass of fluid at a point has a gravity component and a component caused by the gradient of fluid pressure. Thus, 'Em is a measure of the energy per unit mass or the potential of the fluid at any point of the system. The flow takes place from higher to lower energy levels. For horizontal flow this means that the flow takes place from higher to lower pressure, but this is not necessarily the case when the flow is not horizontal because of the gravitational potential. Darcy's law states, in view of the above consideration, that macroscopically the velocity of the water flowing through the porous medium is proportional to the negative gradient of the piezometric head. The macroscopic velocity magnitude q is a bulk velocity smaller than the actual seepage velocity magnitude qse in the pores, to which it is simply related by the porosity E) as q = eqse° Although Darcy's law was first established by experiment, Hubbert (l956)*, and others, have derived it from the general Navier-Stokes equations for viscous flow. Such derivations stem from statistical considerations and ‘The number in the parentheses refers to the year of the publication. References are listed in the Bibliography at the end of the thesis. 9 simplications of the complicated microscopic flow picture. Although they do not contribute to the formulation of a new law, they confirm the earlier belief that Darcy's law is of the nature of a statistical result giving the empi- rical equivalent of the Navier-Stokes equations. The Navier-Stokes equations for two-dimensional flow of incompressible fluid are: au an an ___ ___1 92 A". 2 at + ax + 3y lo ax + ,0 Wu 2V av av = __i 32 .££ 2 _ 9t+' ax+'w3§ P ay"p er 9. One can summarize the statistical averages by the approxi- mate assumptions 2 u v V7 u = - c-—§ and 'v V’= — c«—§ L L where L may be thought of as being a characteristic length, for example, the pore size d of the medium. Assuming 0.! 6 ___33 .. u — 93¢ v — the Navier—Stokes equations become _aaw ala¢2la¢2___1,ag £6319 2x(at) + ax[2(ax) +2(ay) J - / ax+fd2 3x __§_a¢ a 3:192 $2.92 ___-__1£.2 AEBJ’... 3Y(2t) + ay[2(ax) +2(ay)] {a 8Y+fd2ay g. If,11, /0 are constant, and if I. u xha 10 integration of the equations leads to _ g: . IRS—$2 . (31;)2] + 752 -1“; + gy = f(t). (5) This equation reduces to the Darcy law in the case of steady flow provided the inertia terms in the square brackets, here representing the kinetic energy of the fluid per unit mass, can be neglected. Equation (5) with 2.9 at: then reduces to set equal to zero and f(t) constant for steady flow Q = (p + 3'y) + constant .K— J“ which is equivalent to Equation (2), if we identify k with K2’ x41 ° usually presented in a rigorous approach, it is sufficient Although the above derivation omits some arguments for present purposes. For the detailed discussion of these equations, see the book by Muskat (1946). For the problem under consideration here and outlined previously, the governing equations and boundary condition are summarized as follows: The two-dimensional Darcy's equation in an isotropic medium is —-> K C q -:“—(Vp- fgjh (6) The continuity equation in general form for a porous medium is V°F5=—9—ién <7) ’See, for example, Collins (1961). 11 For incompressible fluids, /7= constant. Then Equation (7) becomes V’c'f: o. (8) The boundary conditions relevant to the problem are described in the following: At solid boundaries normal velocities will vanish, which implies that the normal derivative of w must vanish there. Therefore, since w = kh, we have .‘3’ v = O and -§L— = 0 (Bl) n n where n indicates the direction normal to the solid boundaries. If there is a sharp interface between region 1 and region 2, the pressure on both sides of the interface will be the same. Therefore, the following expression for the elevation Y of a point on the interface may be derived by equating the two expressions for p1 and p2 obtained by solving h ='—2 + Y for p. r Y = 72 2:1 71[ :: h2 '- hl] or Y = fi£§f r12 - hl] (132) where F— 3'2 - Y1 12 If the liquid 2 is under hydrostatic condition, the total head h2 will be constant in region 2. In this case, if we take region 2 for the reference state, we may set h2 = 0, and Equation (B2) becomes Y=-—1ho (B3) 1" 1 Also, normal velocities are the same on both sides of the interface, and we have (V ) = (V2)no (B4) C. Review of Literature In the study of flow through porous media, the location and the movement of the interface between two fluids are of particular interest to hydrologists and petroleum engineers. Badon-Ghyben (1888) and Herzberg (1901) stated the hydrostatic equilibrium condition for a fresh water lens floating on top of salt water in a porous aquifer. Since then several researchers have investigated the more relevant case that either one of the two fluids or both are in motion. The most common assumption among researchers in this field is that the salt water is static. Todd (1953) and Kitagawa (1939) combined this assumption with the one-dimensional form of Darcy's law to determine the position of the interface as a function of steady fresh-water discharge. This is similar to the case of gravity seepage, the problem in which the Dupuit (1863)- Forchheimer (1886) theory is utilized. Hubbert (1940) l3 established a general description of the flow of two fluids at either side of a steady interface between fluids. For the case when only one fluid is in motion, the hodo- graph method has been applied to several particular cases. A "hodograph" is a representation of a dynamical system in which the co-ordinates are the velocity components of the particles of the system. Noteworthy of investigators reporting on this method are Henry (1959), Glover (1959), and Kidder (1956). Several major contributions in the analysis of unsteady flow with emphasis on the movement of interface are those of DeWiest (1959), DeJosselin DeJong (1960), Kidder (1956), and Bear and Dagan (1964). DeWiest, in his work, discusses the gravity flow with a free surface. The idea of considering the unsteady flow as a perturbation in time of the final steady state is used. The perturbation velocity potential satisfies Laplace's equation in a dimensionless hodograph plane. The boundary-value problem is then set up mathematically and a numerical example worked out. DeJosselin DeJong based his work on the con- cept of replacing the two different fluids by one hypo— thetical fluid and of treating the different fluid properties by use of singularities along the interface. By knowing the position of the interface and the boundary conditions at a certain moment, the subsequent motion of the interface can be computed. He gives a mathematical example for a two-fluid system contained in a confined l4 infinite aquifer of uniform thickness with interface initially at a vertical position. Kidder treats the case in which the interface motion is caused by gravitational force acting on the two liquids. A general solution is obtained by direct potential theory methods and employs a method of approximation similar to that used in the linear theory of water waves. The method is applicable only when the initial slope and curvature of the interface are not too large. In a very recent publication, Bear and Dagan derive approximate expressions for the movement of the interface in a confined coastal aquifer caused by a sudden change in the rate of seaward flow of fresh water. The solutions are based on the Dupuit assumption. The range of validity of each solution is determined by comparing it with results of experiments. In the present work, the solution to the problem of determining the movement of the interface in the tran- sient state is obtained by successively making use of the one-dimensional flow assumptions and the two-dimensional flow assumptions. In one-dimensional flow, the governing differential equations are replaced by finite—difference equations, and these are solved numerically. In two— dimensional flow, the interface is considered as a dis— tribution of sources along it, and the concept of Green's function is utilized in solving the problem. A case of finite aquifer of rectangular shape with interface ini- tially of the parabolic shape is solved. The two liquids are assumed incompressible and immiscible. CHAPTER II ONE-DIMENSIONAL ANALYSIS A. Derivation of the Equations For the first approximation one—dimensional flow is assumed, since the vertical component of the velocity is negligible in the major portion of the flow field. The simplifications which are made in this part of the analysis are: incompressible flow, sharp inter- face, the specific geometry shown in Figure 4, and an isotropic porous medium. The profile of the porous medium is of rectangular shape with length b and depth d. For the example worked out-g = 20 so that the interface in Figure 4 is much more nearly horizontal than it appears to be in the figure, and the one-dimensional flow assumption therefore seems reasonable. For this dis- cussion y is measured positively downward as shown in Figure 4. Hence, in the equations of Chapter I, y must be replaced by d — Y. The two fluids are of specific weights Y 1 and 72.. I The initial interface position Y0(x) is assumed to be that corresponding to equilibrium conditions with the fluid in reservoir 1 at a level of a0 above the top of the rectangular region, while the level in reservoir 2 15 l6 _._.2_._._. A b I] v a0 a x I ’ [/111/111/‘4111111/111" ___L__ L - ‘ . ~ . ‘ ‘ \r _ ~\ ~ ‘ ‘ r77“ ~fluid l —‘ r / reservoir 1 ’ P ‘ ' reservoir 2 ( : - f‘ -(’ -- d a ‘ ’ f‘ ( :- a —.—r 7' Y - - f g _ - , , oo ,— , ’ ‘ r .‘ . ‘ r . 1 1 - ‘ Y / ’ 4 , 07:, f’ - ’ . , ,fluid 2 ' a ’ , / ’ t . o ‘ ’5’ I / ” I ITIII/I ‘- Ill IIFT77 X=X2 7 I I I 11’] I l I / I // Figure 4.--Flow Geometry for One—dimensional Analysis is zero. Then at t = 0, the level of reservoir 1 is suddenly lowered to the new level "a," which then is maintained while the interface moves toward region 1 (upward and to the left in Figure 4) approaching the new equilibrium position Y°o(x). the transient motion of the interface. The problem is to investigate Because the one- dimensional analysis appeared to be inadequate, only the early part of the transient motion was worked out in the numerical example to be presented below. The applicable equations are Equation (1), which may be written as —+ q "'th, (9) 17 and Equation (8), which is repeated here V“; = o. (8) The one-dimensional flow analysis supposes that in either of the two regions the flow velocity is horizontal and independent of y. If 01 is the total volume flow rate through a cross section x = constant of region 1, con- sidered positive for flow to the right, and u1 is the horizontal velocity component there, positive for flow to the left, then 01 = -ulY since Y is the distance down to the interface, or Q1 = -uld at a section to the left of the point where the interface intersects the bottom of rectangle in Figure 4. Similarly the flow rate and velo- city in region 2 below the interface are related by 02 = "(d - Y)u20 Applying the assumptions indicated above and shown by Figure 4, Equation (9) becomes in the zone 1, and we have 01 = k Y'-——u (10) If zone 2 is assumed to be static, the boundary condition (B3) of Chapter I, part B applies, whence db 1 3’2- 3’1 dY dY ai- = y1 a; = Ta)? (11’ 18 (The sign change is necessary because Y is now positive downward.) When this is substituted into Equation (10), that equation can be integrated to yield 2 20 Y = 1 x (12) I“? if we impose the condition that Y = 0 at x = O. This is the equation of the steady interface for any specified value of 01' Substituting Equation (B3), with a sign change, into Equation (10) and integrating with respect to x, we have Q — .l.El.E£f for 0 é x e x l " 2 I“ x 0 where we have used the condition h = 0 at x = 0. For 1 x 7>x0, integration of Equation (10) for Y = d = constant with x = b when h1 = a, gives Ql(x - b) = kld(h - a) or k1d(h1 - a) = ‘é (é Ql x _ b for X0 x b. If these two expressions for Q are equated with x = x 1 1’ noting from (BB) that h = I‘d, the following relation 1 between the level "a" in reservoir 1 and the co—ordinate x , where the transient interface intersects the bottom of l the aquifer, is established == 3%(1 + JR). (13) .2 d xl 19 The resulting Q1 for this relation is, after solving for ‘EE and substituting into one of the equations for Q1 1 _. .1. 9. .19.. 01 " 21“1‘1(1:)17xl — k (3m - 3k (33) (14) 7 l d 2 l b ° Darcy's equations are: l9h1 Q1 - le-b—JE- (15) 9% 02 = k2(d—Y) ax . (16) The continuity equation for region 1 is, after integrating two-dimensional form of Equation (8) from the upper boundary to the interface, Y ggul £9v1 a ‘Y Jo(—a X + a y)dy —- - efijody = temporal rate of decrease of fluid caused by inter— face movement which then reduces to u (iiY)+u13:+9%—%=o (17) where the bar indicates a mean value. u For one-dimensional flow, 7§_% = constant across the cross section and since 01 = Yul ‘See p. 35 of "Advanced Mechanics of Fluids" by Hunter Rouse (editor). 901 aY a x = - a -a—E o (18) Similarly, for region 2 302 BY 3 x = + 6 E o (19) Adding (18) and (19) and integrating with respect to x yields + Q = f(t). Q 2 1 Alternatively, a definite integration with respect to x from 0 to b yields, after substituting (14) through (16) x1 hl h2 b h J01 [1(1ng + k2 (d-Y)aax2 de +J k (ii—3'- dx x1 1 39x d2 1k .— 1 + 02 = 3k l-EJ7[I‘d — l] — constant. (20) Equation (20) indicates that QT remains constant during the transient state. The expression in Equation (20) is equivalent to that in (14) as can be seen by solving (13) x for'—% and substituting in (14). This is to be expected, since at xl'4 x é’b the entire discharge is in zone 1. B. Numerical Example As a specific example, we will consider the porous medium in Figure 4 with the following properties: X l§-=2o, —-11-)=o.5, ..%=o°95, I‘= )2—2’131_ (71 4O ' 21 Ky Kar We takedxl = ° then since k = 1 and k = 2 1 x42’ ’ ' 1. X11 2 A2 we have k2 =-§;gkl. Initially, an equilibrium condition 1 a prevails and Equation (13) yields-—% = 1.517 corresponding to x1 = 0.5b. The elevation a0 is then suddenly lowered to‘% = 1.02617 corresponding to an equilibrium interface with x2 = 0.95b. At the initial steady state the flow _ l . rate is QO --§6kldI" according to Equation (20), and during the transient state and at the final steady state _ 1 QT " 38kldI-‘ 0 Introduction of the following dimensionless variables allows for considerable simplification: Q Q 1 2 r . x 521%., 512%", Y fit " =3, T T h’— klhl hi— k2h2 t, — QTt - , -- , "" O 1 QT 2 QT 9d2 Equations (15), (16), (l8), (19), (20), and (B2) become: . ah' .f21 = y' B E' (21) X h/ 322 = (1-y’)9,2 (22) 93c 399”} = —Y—: (23) x t 352,2 = 7—3-17 (24) x t 511 + 9.2 = 1 (25) 22 ’-QT (' ’) (2) Y fi—qu—‘hl-hzo 6 The initial and final positions of the interface are computed from Equation (12), using the values of Q0 and QT respectively, given above in the paragraph preceding Equation (21), as: I2 ’ Y0 = 0.1x (27) y’ 2 = o oszex’ (28) w o 0 Immediately after the reservoir surface of liquid 1 is changed, the value of Q1 on the interval b a X'a 0.5b is the final Ql = QT because Q2 = 0 here and Q + Q 1 2"QT“ constant. Hence.521 = l on this interval, and since y' = l on this interval, Equation (21) integrates to yield (hf)0 = x' + constant. Therefore, we have for the initial distribution of hi to the left of x = 0.5b immediately after the lowering of the liquid 1 level to its final value: (h ) = x’ + 19 0.5b é-x'S b. (29) I The distribution of hi and h2 to the right of x = 0.5b is obtained by substituting (21) and (22) into (25) together with (27) yielding, after substituting the numerical data, I 2069 o I ‘ 4 (h ) =-———-./§' - 0.0455x o - x - 0.5b (30) 1 o {35 and (hé) -0.0455x’ o E x 5 0.5b. (31) 0 ll 23 Computed results of (29), (30), and (31) are plotted in Figure 5. C. Finite—difference Solution Procedure The difference equations corresponding to Equations (21) through (26) are: (h’) - (h/)° (£11)1 j = y; j 1 i+l,j , 1 l’laj (21A) 3 ’ 2(AX ) (h’) - (h')‘ (522)1 j - (1 - Y; j) 2 1+1’j 2 i-l’j (22A) : ’ 2(AXI) (521)i+1,j ' (522’1-1.j = yiaj+1 ' YiAJ (23A) 2m)“ At’ 2(Axl) At, (S)1>i,j + (522)i.j z 1 (25A) , QT ,. ' Ym’ . gammy“ ““ ”12):.9" (26A) In these equations, the index i refers to the x' increments while j refers to the t' increments. The numerical computations are carried as follows: 1) Choose suitable intervals for Ax’ and At’. 2) From Equation (23A) compute, starting with j = 0 24 , _ (521)1+1,j“ (£21)i-l,j +y: ,j+1 _ I 13:) o 2(9E7) At 3) Substitute Equations (21A) and (22A) into Equation (25A), replacing j by j + 1. Equation (26) is then substituted into the resulting I equation to eliminate (h1)i+l,j+l and (hl)i-1,j+l’ yielding I / (h2)i+l,j+1 - (h2)i-1,j+1 = 1 n 0095 ’ o 2(Ax') Ax’ 1’j+1 I I (Y1+1,j+1 ‘ Y1~1,j+1)° 4) The equation obtained in step (3) is then substituted into Equation (22A), giving (311)1,j+1 = (yi,j+l) ° 0.95 [1 *‘Z;" (yi+l,j+l ” Yi-l,j+l) (1 ' Y1,j+1)] ° 5) (I? ) is then obtained from Equation (25A). 2 i,j+l 6) Repeat the steps(2) through (5) for new increments of j. D. Results and Discussion A numerical example with Ax' = l and At, = 1 was carried out, and the resulting interface locations corres- ponding to t' = 1 and t' = 2 are shown in Figure 6. This 25 increment in x’ is equivalent to dividing the entire length of the rectangle into twenty intervals. Since the initial interface occupied only half the length, it amounts to dividing the initial interface into ten equal horizontal intervals. The time interval was selected arbitrarily, to be checked for reasonableness by requiring that the calcu- lated motion of the interface during the time interval be small compared to the distance from the initial equilibrium position Y0(x) to the final equilibrium position Yoo(x). These computations are designated as Case 1. Observation of Figure 6 reveals that the computed values of Y tend to go above the final steady state near the right-hand end of the flow field. This is not likely to be the case, since oscillation is not anticipated. The difficulty may be attributed to the relatively steep slope of the interface in this region, where the vertical component of the velo- city is not negligible and therefore the one-dimensional analysis is not a good approximation. An attempt to improve this discrepancy was made by taking smaller inter- vals (Ax’ = 0.02, At, = 0.02) in the difference equations. The results of this as shown by Figure 7 disclose a similar tendency. Here, the locations of the interface computed for t = 0.02 and t = 0.04 are plotted. This is designated as Case 2. Another calculation was made after interchanging the initial and final positions of the interface, i.e., making the final steady-state position of the interface in the previous case the initial position of 26 the interface in this case. The result of this for t' = 1 and t’ = 2 as given by Figure 8 indicates that the values of Y tend to fall far below the final steady state near the right end. This is designated as Case 3. Again, we are unable to overcome the difficulty mentioned in the previous case. All the data pertaining to Figures 5 through 8 are listed in Tables 1, 2, and 3. The results obtained in these analyses lead to the conclusion that the one-dimensional approximation cannot be valid near the outflow end, though it may possibly be a good approximation elsewhere. The two-dimensional analysis will be used in the next chapter to yield a refined solution. s\x .m> a\n mo Loam m .hSHAW .o\x lllv o m .4 ,6 .w 8 m. 1: 8 m: 8 ml..- 3, 3:.13 - 1.- . _- :TQIiL--.!.-:-.|er.-.}_?--!. . Ll 4_ q A L m . 1v _ \.\ U (\ III.“ _ mmmnuumomMMmch Mo pcnmm>ou U\N Ill! OF 0 P (\1 - . .. @ madmam 3- 3 ) 4. N H .le.l|l.-|ll. —H .9. (lift @cmmmq m ammouumommpmch Mo phmumboa b madmam N 0 mo. 0?. dm. mm. Cd. 0. o 2 .-Itqu- I . I - T. . : .-.I-..--4.:.I -. . :. ‘1 I I I I I , \Innx .\ I I y 1 \ ..HMMHHL\\ (\...\(.MMI : 3 - 1 - . (\fiflvkfl\\ ,IIIIII Awaw . 40.0“” w.I I.Iunu \.\\ \\(\\.(\\\\I I)‘ I a . T \ III. I 2 . \n\ X . No.0 I .p IIIIII \ \ s I. ..I _ \\\WVE\ unomoq ..\ \\\ \ " \\\\\,\ u .\\\ R O I (III: II II III I I-.. I II .............. I I IIIIII 3.0 \O m mmmonlmOMMpmpmH mo pmmam>oz 31 m u u now coauwmoa mummumuca 0:» ma N» H u p now cowuflmoa mummumucfi m£# ma a» mwoz H¢.w Hq.m om mm.@ mm.@ ooo.H mH‘ oo.m mo.m mum.o mH om.m mm.m mHm.o 5H v>.m mb.m me.o mH mm.m mm.m mmm.o mH mm.m mv.m mmm.o ¢H mH.m mm.m bmm.o mH mm.¢ OH.m vmboo NH mb.« wo.¢ 005.0 HH bum.o mm.v 00¢.Hu hb.¢ mmn.o ooo.H OH oom.o om.v mHm.Hu Hm.¢ mmm.o mvm.o m mom.o vvm.o mo.¢ OFH.HI mm.¢ m¢m.o mmm.o m omu.o mmh.o om.m mNo.Hu mH.¢ mom.o 5mm.o p mmm.o me.o m¢.m mfim.ou mm.m mvm.o mbb.o m Hum.o Nvm.o om.m cm».on mm.m mHm.o For.o m Mbv.o 0mm.o hm.m mam.ou mm.m mm¢.o mmm.o H mmm.o vmv.o mv.m wm¢.ou mm.m 5mm.o mHm.o m HHb.o Nvm.o mo.m Nam.ou om.m Hmm.o bw¢.o m NHH.o m¢.H mvH.ou mm.H mmm.o mHm.o H o o o o o o o N» HH co m o m o m » ow o\x OH x H av OH x H av OH x H av m m w 924 m mmMDUHm mom ZOHB¢ADm‘yi. The method of solution applies at any time of the transient motion to calculate the next step of the motion when the two fluids are in a non—equilibrium state. In the example to be presented the initial interface position is a parabolic approximation to the known equili- brium configuration for a certain total discharge rate Q. Transient motion of the interface is then examined under the assumption that at t to a higher value and maintained constant. 0 the discharge rate is changed Such a change in discharge rate could be produced physically by raising the level of the fluid in reservoir 1, and this would cause the interface to move to the right. The boundary 36 conditions indicated are obtained from the fact that the reservoirs are under static conditions and that the sur- faces AF and CDE are impervious. Under static conditions the total head is independent of y since the pressure head increases with depth by an amount equal to the decrease in elevation. Thus on FE we have $1 = kla1 and on EC $2 = k Since the normal derivative must be zero at the £232. impervious top and bottom boundaries, we have (<01)y = 0 on ‘ AF and ED and (€92)y = O on DC, where the subscript denotes the partial derivative with respect to y. The boundary condition for w on AB is derived below as Equation (38). l Darcy's law for homogeneous two-dimensional flow through porous media is valid everywhere in the porous medium except along the interface. Therefore, we have h 9h u -k a 1’ V = “k 1, 1 1 ax 1 1 a y (32) 9% ah2 u.2 = “kg—9‘? V2 = “’92-'57 where the subscripts refer to the appropriate zones. If we introduce the velocity potentials $1 = klhl and ¢ = k2h2 into Equations (32) we have 2 gml 9‘91 111 = - , v = ___, 9x 1 BY (33) _ aq’2 9‘92 “2 " " ax’ V2 ‘ ““57 or in vector form 3’ = -VCP. (34) 37 Since the continuity equation must be satisfied, we have, from Equation (8) for incompressible flow, the Laplace's equation 2 K] o = o. (35) Equation (35) applies in regions 1 and 2, respectively, but does not apply across the interface. The boundary conditions along the seepage face AB and the interface BD are derived as follows. The pressure is continuous across these faces. Therefore, the condition q’1 “’2 p=(Tc-;-y)X1=(TCE-y))(2 (36) has to be imposed on them. For the vertical seepage face AB, where $2 = k2a2, Equation (36) becomes )‘2 3’2 1 1 71 2 )1 1 To further simplify the analysis, let /u2 =./Ul°. Then, since K‘X KL?’ 1 2 k = and k = 1 ,1X1 2 ,111 we have‘jigk = k ; and, if we introduce k' = (Jig - l)k , ‘31 1 2 y]. 1 Equation (37) reduces to ‘If the two liquids are fresh water and salt water, respectively, the viscosities are practically the same a 38 w = k - k’y. (38) 1 232 Similarly, for the interface BD Y2 ¢ = ¢ - k 0-- 3‘1 1 2 1 - l)y = $2 — k y. (39) Also, normal components of the velocities along the interface are the same on both sides and therefore g3¢1 a“) 9n1 9% N (40) since positive n1 and n2 are chosen outward from the respective regions and are therefore in opposite directions on the interface. We introduce the stream function 1p defined by uz—a-Eé v=-aw (41) ay’ 2x which with ¢ satisfies the Cauchy-Riemann equations 23!. ... -22 2.2. 2 .252. (42) 2y ax’ ax 2y° Differentiating Equation (39) with respect to s, the distance along the interface, yields 2Q’1 9¢2 ,9 25 = ‘3‘? -k3—§ (43) where s is defined positive from D toward B. In terms of EP, Equation (43) becomes _Qfi 2 _EE -1621 (44) 89n2 29n2 3 S where we have used Equation (42) with local x and y axes 39 in the s and n directions. Since 2W = — 23L, this 2 9n2 anl becomes 210]. = 2&2 _klaio (45) 3n1 anl as Also, Equation (40) becomes a 1 2W2 l = W . (46) 255 In summary, the problem, in terms of 1;) , is to solve the governing equation V 22f = 0 with the boundary conditions shown in Figure 10. Subscripts x and n denote F WZQ A Region 1 Figure lO.--Boundary Conditions in Terms of 2// partial derivatives. The boundary conditions on Walong EF and AC are obtained from those on w shown in Figure 9, by using Equations (42). These equations also imply that IA 40 z? is constant along the top and bottom. If we arbitrarily set '4’: 0 along the bottom, then along the top IP is equal to the total discharge Q, since along EF-{%¥L is equal to u. The interface condition along DB was derived above as Equation (45). To simplify the problem one step further, a function ‘g7is introduced, defined by 1}): W'% . Then v27? = v 2270 and the problem is to find a function I? which satisfies Laplace's equation in region 1 and in region 2 with the boundary conditions shown in Figure 11. F (#0 A Region 1 Figure ll.--Boundary Conditions in Terms of I? B. Solution Procedure with Green's Function Since a discontinuity exists along the interfaces DB and BA, the governing differential equation does not apply across them. The interface, with normal velocity discontinuity as shown in Figure 11, can be considered as 41 G=0 Figure l2.--Boundary Conditions for G a distribution of sources with sources intensity IP( I, 7) -232- Here per unit length, equal to the discontinuity in 377. :§(s) and 77(5) are the rectangular co-ordinates of the source point specified by the arc—length parameter s. The details of solution for '4)are presented in Appendix V. The results are quoted below in Equations (47) to (49). The solution for If(x,y) in either region is given in terms of a Green's function G(x,y; §,I2) as’ B 1?(xq) = -];4 00 4 G(an; 9 ) = "' Z Z d o E 7 nzbd m=1 n=0 n (48) C05 BLIP-C sin m COS m sin M b d b d 23 + 12. b2 d2 where 04 = 1/2 if n = o d = 1 n ’ n 1 + SnO or 6% n = 1, if n f 0 6n0 = Kronecker delta. When we introduce the source strength f9: kairfi on DB and f): -k’ on BA, Equation (47) takes the form 2} (x,y) B - [D G[x,y; f(s).7Z(s)] k'—S% ds - A k’fB G(x,y; b,7 )d7 II A "k/fD G[x,y;§(7),7]d7 (49) 43 since on DB Y(s) = .7(s). Substituting Equation (48) into (49) then yields I oo nTrx . mny 4k 03 cos b Sin d . ? (X,Y) = - 2 Z >_ 0( n 2 2 n bd m=l n=0 -2_ + m b2 E? d (50) Ocos rm f) )sin-I-nflc—ilz-d7. ‘ If (50) is made dimensionless, we have 1} I: .‘P = - Z Z ”(n COS2 b :1: 2d . 4k d b m=l n=0 n + m 0—) 2 (a) d n (51) 1 [0 cos 21b:- sin Bial- d(-%). J The function 97 in Equation (51) is computed numerically, and the movement of the interface is derived from it. A detailed account of the numerical solution is explained below. The Control Data 3600 Computer of the Michigan State University Computer Center was used for all the numerical computation. The computing procedures are as follows: 1) Assume an initial interface position DB. 2) The interface position is expressed in the form of a group of discrete points (xi, yi). The intermediate points among these discrete points are computed by interpolation using Newton's 3) 4) 5) 6) 7) 44 divided difference formula. This formula is listed in Appendix II. These values are stored in the memory of the computer for later use. The integral in Equation (51) is evaluated along the interface using Simpson's rule with intervals of-égL : 0.001. The results of part (2) are utilized in obtaining the values of-J%. 4) in Equation (51) is then computed along the interface at the intervals of 9% = 0.1 except the points near the end B where intervals of'éfi = 0.01 are used. The numerical data giving 1? at points on the interface is then differentiated with respect to E. d velocity of particles on the interface. dimensionless arc-length to obtain the normal The results of (5) are then multiplied by the dimensionless time increment Ask to obtain the dimensionless displacement 9% of the interface in the normal direction. The horizontal and vertical components of this displacement are obtained by multiplying-Ag by the direction cosines of the normal. The new position of the interface is then obtained by adding these horizontal and vertical components of the normal displacement to the corresponding initial co-ordinates of the interface. (Note that these x and y components are not actual x and y 45 components of the displacement of material particles, since particles would also have some velocity components tangential to the interface, computed from the normal derivative of E/ and therefore discontinuous at the interface.) 8) These co-ordinates are then substituted again in Newton's divided difference interpolation formula to obtain the intermediate points. 9) Parts (2) through (8) are then repeated to obtain the subsequent positions of the interface. In Appendix III, the entire computational schemes are presented in the form of Flow Diagrams so that the computer instructions can be easily programmed from them. Two sample FORTRAN Programs are presented in Appendix IV for reference. C. Numerical Example For the purpose of illustration, a numerical example, with the following data is presented: b a» = 20, —'Q— - 300 The initial interface position was obtained from the equilibrium position calculated by Henry (1959) for a discharge rate represented by-Eég : 40.025. Co-ordinates of eleven points along the initial interface are shown in Table 5. Additional points on the interface were obtained by interpolation as described in step (2) of the procedure given in the preceding section. At time t = 0 it is 46 assumed that the discharge rate is suddenly increased and maintained during the transient motion of the inter- I face at a rate represented by-Eag = 30. D. Results and Discussion Table 4 shows the values of 9? (labeled PSI) between-é = 0.1 and3§ = 0.8 along the initial interface, the results of step (4) of the first computation cycle. The number k in the table indicates the number of terms in the double series of Equation (51); thus k = 60 means a partial sum with m going from 1 to 60 while n goes from zero to 60, a total of 3,660 terms in the double sum. Since the results for k = 50 agree with those for k = 60 to approximately four significant figures, it is hoped that sufficient convergence was obtained to give results accurate to approximately four significant figures. This result is plotted in Figure 13, showing that the relation between. 'LIJIand-x is linear except at the points beyond d «é = 0.80, which are computed separately and given in Table 6. Two cycles of the calculation outlined in part B were performed to obtain the interface positions at times corresponding to-EEL equal to ten and twenty. The resulting interface positions are plotted along with the initial interface position in Figure 14, based on Table 5, Which gives the position co-ordinates of eleven points on the interface after the motion. It appears from Figure 14 that the motion during the second time increment was 47 considerably smaller than during the first. This plot shows a more reasonable motion of the interface than was obtained by the one-dimensional analysis. However, a closer examination of the region near the right-hand end of the flow field shows still some discrepancy. In the region between‘fi = 0.8 and«% = l, the computations were performed at intervals‘éfi = 0.01. The results for E? are given in Table 6 with the series carried out to k = 110 in order to get sufficient convergence. The results are also plotted in Figure 15, and it is seen that the linear plot of Figure 13 extends to about-x = 0.89. Beyond this point the results become d erratic and unreasonable. The position of the part of the interface between-é = 0.8 and-é = 1 after-Eg— = 10 is shown in Figure 16, based on Table 7. Near the right- hand end the large displacement shown is in the region where the solution is erratic. Perhaps a better result could be obtained by taking a smaller interval in the interpolation formula and the numerical integration. Except for this erratic behavior very near the end, the two—dimensional analysis predictions seem reasonable, and the unreasonable behavior occurs in a much smaller part of the field than was the case with the one-dimensional analysis. ‘48 ..\» 53% mo 3: 2 933.. Allin: ox» m.c '96.. « — mo.1 3.... no.1 Nc.n ado: ‘ 1 I II!!! II. I All. Ill-Ill" llllllllllllllll [IL-I ON é omuotuooduuoan no anoaobo: a— ousmum AIU\N ..N— N— O-. 49 at!) \0 Hr q H _ 4 H . o_ u a\.xp Hove» ooauuoan HaHHHuH "IN 0 xxx. cannon N. (T‘) \ f) 1.1 Y/d—* 9 “I . \J \C) (I) .34 (\I U ) r) 50 .N"\ ’3 [l N" ( ) "\ -" r- . 1 0 v .4 51, c.o~ ”.mn m omaunuounmuoucu we acolo>oz 0H ousuum 3x oH . a\.xu page. aawuwcw vsouoq o.~ e\» TABLE 4 K320 VALUES OF PSI ALONG INITIAL INTERFACE FOR Y/DIOoI T0 0.8 Y/D 0010 0.20 0030 0040 0050 0060 0070 0080 0.10 0.20 0030 0040 0050 0060 0070 0080 0010 0020 0030 0.40 0050 0060 0070 0080 0010 0.20 0.30 0040 0050 0060 0.70 0.80 0010 0020 0030 0040 0050 0060 0070 0.80 52 (TABUEATION FOR FIGURE 33) PSI ~0553055545737E-02 “0115720761094E*01 ~0176355510117E*01 ~02364594646575-01 -.2955960623OOE‘01 “03511790625855‘01 ~0398317114792E‘01 -0448916755369E“OI -.582795472379E’O2 “.119473692097E-01 “0180447219584E‘01 ~0241584277297E-01 ~03018955400555‘01 ~036099911249OE‘01 “04140437485775-01 -04536879646625-01 *05940177029125-02 ~01206897914O7E-01 -0181970215872E-OI -o2431730682155~01 ~030396296819SE-01 -.363412571118E~01 -o4l96951663IBE-OI “04712934398195-01 -0599469170993E-02 ~0121334O403léE-OI ~01827048674865-01 ~02440304063355-OI -.305048160917E-01 ~0365143438139E—01 -o422962561803E-01 -o474693501048E-OI ~06028487584IZE-02 -o121729893888E-OI -o18316377854OE-01 ~0244519188225E-01 ~0305669839261E‘OI ~0366O901052355-OI ‘04247035935118-01 -0477747216275E-01 DIFFERENCE -0553055545737E‘02 *0115720761094E~01 “0176355510117E-01 -.236459464657E-01 -.29559606ZBOOE-01 -.351179062585E-01 -0398317114792E-01 -0448916755369E-Ol -.297369296528E-03 ~0375293100216E*03 ~0409I7OQ46587E’03 -0512481264217E-03 -.629947775480E~03 -o9820049909IBE-03 -0157266337832E~02 -o4771209287SOE-03 -01122223053SZE-03 -0121609931060E*03 -0152299629OISE‘03 ~0158879091494E*03 -.206742814506E-03 “0241345862385E*03 -.5651417741615‘03 -0176054751595E—02 ~0545I46807563E-04 -.644248907S4aE-04 -o734651612234E-04 ~08573381228635~04 ~0108519272093E‘03 -ol730867024965~03 --326739548669E*03 ~0340006123239E-O3 -.33795874I744E~04 -.395853574OOZE-04 -o¢58911054018E-04 -o48878188864ZE-04 -o621678341326E-O4 -o9466670963015-O4 ~0174IO317IOSSE-03 -o3053715227I4E-03 53 TABLE 4--CONTINUED Y/D PSI DIFFERENCE K860 0010 ~0605041735282E‘02 -.219297686504E-O4 0.20 ~0121992332453E-01 -.262438561556E-04 0.30 #0183459515629E"OI -o29573708617OE‘04 0.40 ~¢244866174082E-OI -o346985852964E-O4 0.50 ~0306114952320E“01 -o445113055328E-04 0.60 -0366647209405E-01 -.557I04167470E-O4 0.70 ~0425997797727E“OI ~0129420421224E-O3 0.80 ~0481022714914E‘OI -o327549862974E-O3 NOTE K INDICATES THE NUMBER OF TERMS IN EQUATION (51). FOR EXAMPLE. K860 MEANS THAT 60X60=3600 TERMS ARE CARRIED OUT IN THE SUMMATI! TABLE 5 AT T=0 AFTER T=10 NOTE ... 1. 2. 54 MOVEMENT OF INTERFACE POSITION--CASE 4 (TABULATION FOR FIGURE 14) Y/D 0.000 0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900 0.982 .19738456E .29700272E 039645103E .49560301E .59399079E .70233221E 081574804E (AT T -.22710731E-02 .97638290E- 01 00 00 00 00 00 OO 00 THE DATA FOR SECOND SET OF VALUES CORRESPOND TO THE FIRST 9 VALUES IN THE FIRST SET CONTROL DATA 3600 WAS USED FOR T = IBM 1620 WAS USED FOR T = 10 X/D 0.000 .380237352438E+01 .720449719549E+01 .102063710230E+02 .128079950071E+02 .150093691478E+02 .168104934459E+02 .182113678999E+02 .192119925101E+02 .198123672759E+02 .200000000000E+02 .61477886E-O4 .36941965E 01 .71083374E 01 .10122245E 02 .12735916E 02 .14949369E 02 .16762655E 02 .18174974E 02 .19186239E 02 (AFTER T= 0) 0 AND 10) TABLE 6 VALUES OF PSI FOR INITIAL INTERFACE POSITION Y=O.BIO T0 0.982 K8 20 Y/D .810 .820 .830 .840 .850 .860 0870 0880 I 0890 0900 0910 0920 0930 .940 .950 0960 0970 0980 0932 0810 0820 .830 .840 0850 0860 .870 0880 .890 0900 .910 0920 0930 .940 .950 .960 .970 .980 .982 SS (TABULATION OF FIGURE 15) PSI *0452184037582E-OI “.453708840054E-01 ‘0453144403000E‘OI “.45017956219IE‘01 -.444540875447E“OI “o435993889711E‘01 ‘0424344879320E'OI ‘o409443807417E‘OI -0391188634036E‘01 ‘0369530592252E‘01 -.344483332432E‘OI. -.316116333968E‘01 -.284569470103E‘OI -.250050968927E“01 “.212836335268E“01 “.173264630052E‘01 -.I31732188562E-01 -.886840445251E“02 ~0799351987516E‘02 “0457766126623E‘01 *0463284042831E“OI -0469820578153E*OI ~o476526464969E‘OI ‘0482256676085E‘01 ~.485752950466E*01 ‘0485827919299E‘OI “.481503599367E‘01 -.47207453694IE‘01 -.457093757883E-01 ‘043632142213OE“01 -o40960079123OE~01 -0376849203712E*OI ~.3380444I7331E“01 ~.293274613025E‘01 -.242819787385E-OI -.187233390643E-OI ~.127392746388E‘Ol ~0115018000693E-01 DIFFERENCE -.452184037582E-OI -.4537OBB40054E-OI -.453144403000E*OI -.450I79562191E-01 -.444540875447E-OI -.4359938897IIE~01 -.4243448793205‘Ol -.4094438074I7E-OI -.391188634036E‘01 -.369530592252E~OI -.344483332432E-Ol -.316116333968E*OI -.284569470103E901 -.250050968927E‘01 -.212836335268E-OI -.I73264630052E-OI -.I3I732188562E~01 -.886840445251E-02 -.799351987516E‘02 -.5582089043025-03 -.95752027756OE-03 -.166761751556E-02 -.263469027785E‘02 -.377158006362E*02 ~0497590607SBSE‘02 -.6I4830399820E-02 -.7205979I9470E~02 -.808859029086E-02 -.875631656367E-02 -.918380896939E“02 -.934844572679E-02 -.922797336098E-02 -.B79934483979E-02 -.BO4382777569E-02 -.695551573343E-O2 -.555012020792E~02 -.387087018644E-02 -.350828019422E-02 Kg 30 K: 40 Y/D .810 .820 .330 .340 .350 .860 .870 .380 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 .810 .320 .830 .840 .350 .360 .870 .880 .890 .900 .910 .920 .930 .900 .950 .960 .970 .980 .982 56 TABLE éfi‘CONTINUED PSI ‘.472920946682E‘OI ~.4735233B6812E~01 “.474514539543E-01 -.477125641155E‘OI -.48186144323OE‘01 -.488250453061E‘OI -.494889206624E‘OI “.499733504439E‘OI -.500545813906E‘01 -.495342207585E“OI -.482702450587E*OI -.461712520708E~01 ”.431911012289E*01 ‘.393110517856E‘Ol. ‘.345345048154E‘OI -.288942944855E‘01 “.224660673521E‘01 -.153785767827E*Ol -.138980940221E‘OI -.4aoa713677ose~01 -.4as735333.595-01 ~.48803333819SE-01 -.43319..14735£~01 -..aao970965105~01 -.48977404eoooe-ox «.4942710758585-01 -.500631167379E-01 ~.5069482850528‘01 -.5098146123645~01 ~.5057925263345-01 ~.492239518324E~01 -.4678353433855‘01 ~.432035788072E~01 -.3644187949.3£-01 -.324799319184£-01 ~.2539810023308-01 -.174167838861E-o1 -.157412018762E-01 DIFFERENCE -.I51548200578E-02 -.10239343983IE-02 -.469396139124E‘O3 -.599176182746E-O4 .395232855233E-O4 -.249750259176E-03 -.906128732822E-O3 -.182299050724E‘02 -.284712769651E-02 -.382484497008E-02 -.463810284607E-02 ~.521117294746E-02 -.550618085806E-02 -.550661005327E-02 -.520704351235E-02 -.461231574678E’02 -.374272828783E-02 -.263930214380E-02 -.23962939527OE-02 *.795042102254E-03 .-.122119466416E-02 -.I351879865IIE*02 -.110687735834E~02 -.623565328176E-03 -.152359293683E-03 .7181307683065‘04 -.897662939693E~O4 ~.640247II44685-03 -.l44724047826E*02 -.230900762475E‘02 -.305269976168E‘02 -.3592433I0966E-02 -.389252702102E*02 -.39073746794OE—02 -.35856374325SE*02 -.293203288136E-02 -.2038207IO343E~02 -.IB43IO785382E-02 K2 50 K8 60 Y/D .810 .820 .830 .840 .850 .860 .870 .880 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 .81 0 .820 .830 .800 .850 .860 .870 .880 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 57 TABLE 6--CONTINUED PSI -.481634665775E-01 ‘.487I47646396E‘01 *.493299794303E‘01 -.497661655689E“01 -.498873843368E‘01 ~.497962164540E~01 “.497808570639E‘01 -.500798562876E“OI ”.506720480247E*OI ‘.512978024577E~01 “.515718601426E“OI ‘.510197363445E*OI ‘.492309120069E-OI -.460339127058E‘OI -.4l4321020340E-01 -.353883813827E-OI *.278702420605E‘01 “.191220144316E‘01 -.172729686939E‘OI -.485827322533E‘01 ~.489249962920E‘01 *.4935605986685-01 -.499627180303E-01 “.505022735364E‘01 *.506885582625E-01 “.505629894935E-01 -.504551362523E~OI -.50692827304OE-01 -.512968505820E*01 -.518915871333E-01 -.519284426817E‘01 ‘.508612640420E‘01 *.482110I96281E‘01 -.43817778017SE-01 -.377589277676E-01 -.299859439787E-01 -.206186886688E‘OI -.1861471868I4E‘01 DIFFERENCE '-.7632980669005~04 -.14123129404IE‘03 -.526645610997E“O3 -.946724095236E-03 “.1077674685825-02 9 . 8 I 88 I I 854522E“ 03 -.3637494783185*03 -.167395501197E-04 .2278048123155’04 -.3I63412211585*03 -.9926074590183‘03 -.1795784512115‘02 ‘.244737766847E“02 -.2830333898675‘02 -.299022253974E‘02 “.2908449464525‘02 -.247214182724E‘02 -.I70523054523E‘02 ‘.15317668I8088‘02 -.4I9265676101E“03 -.210231652090E-O3 ~.2608043632865-O4 -.I9655246I603E~03 -.6I4889199831E-03 -.892341808139E-03 -.782132428867E-03 -.375279964835E-03 -.207792791120E-O4 .951875335901E-06 -.3I972699OTIBE-03 -.908706337214E-O3 -.163035203517E-02 -.ZI77IO692235E-02 -.238567598357E-02 -.237054638459E-02 -.ZIIS701918265*02 -.I49667423742E-02 -.134I74998722E-02 K: 70 K: 80 Y/D .310 .820 .830 .800 .850 .860 .870 .880 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 .810 .820 .830 .840 .850 .860 .870 .880 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 58 TABLE 6--CONTINUED PSI “.487233091830E‘01 “.492770942314E‘01 “.4963586672095‘01 ~.500081645529E-OI ".505883439895E‘OI “.511448045057E‘01 -.513168679303E‘OI ~.511243153218E‘01 -.509902017104E-01 *.512766537460E‘01 -.519274954102E-01 -.523762289595E‘OI “.518802227132E-01 ‘.498314128214E*OI “.458036090546E“OI -.397633889501E“OI “.317970749391E‘OI *.219¢85890513E“OI ‘.l98128739656E“OI ‘.488057376664E‘01 -.493395512865E“OI *o499l60038438E‘OI “.502793829161E“OI -.506330617296E‘01 -.512172595030E“OI -.517383527884E“01 -.5179I6155295E‘01 ~.515100689212E‘01 -.514382958980E“OI -.518859315431E‘01 -.525279511159E‘01 -.525285820899E‘01 -.510266760015E-OI -.474280648501E‘01 ‘.4ISI39840159E‘01 “.3338314104485‘01 “.231252241239E‘01 “.208794458089E‘OI DIFFERENCE ".I405769298845‘03 “.352097939875E“03 -.2798068544588‘03 ’.454465225637E‘04 -.8607045219985“04 -.4562462436305‘03' '.753878436633E‘03 -.6691790695335-03 -.297374407644E*03 .2019683688605-04 -.359082769144E~04 ’.447786277618E‘03 -.1018958671035‘02 -.1620393193705‘02 “.1985831037055‘02 ”.2004461182915‘02 “.181I130960098‘02 ‘.I32990038266E‘02 -.1198155284358‘02 ‘.824284834394E~04 ~.6245705¢81125~04 -.2801371229075‘03 -.2712183632035‘03 -.447I77408212E-04 -.724549963721E-O4 -.421484858634E*03 -.667300207I77E-03 -.519867210649E-03 -.161642152310E-03 .415638678540E-04 -.I$I722156264E-03 -.648359377024E-03 -.II9526317942E-02 -.162445579557E‘02 -.1750595065505-02 -.158606610607E-02 -.II7663507264E~02 -.106657184326E-02 K: 90 K3100 Y/D .310 .320 .830 .840 .850 .860 .870 . 880 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 .810 .820 .830 .800 .850 .860 .870 .880 .890 .900 .910 .920 .930 .940 .950 .960 .970 .980 .982 59 TABLE 6--CONTINUED PSI -.490242524080E‘Ol “.4943264OIIO9E~OI ".499569622698E-01 “.505279814039E-OI “.508639785854E‘OI ‘.512372014I71E‘OI *.518449229319E‘OI -.52245044567IE'OI ‘.520947900208E*OI ~.517753752758E‘OI -.519113687056E-OI -.525386022244E‘01 -.528972172862E-OI “.519190544233E‘01 “.487641899512E-01 -.430385188425E“01 -.348078925490E“OI *.241740687696E-OI “.218306809391E‘01 -.49o7034238425-01 -.496302525367E-01 -.sooz97199484£-01 -.5057690808468~01 -.511102328077E~01 -.514012544823£-01 “.518433598423E‘OI -.524399327172E‘01 -.525873816587E‘01 -.522253125819E‘OI -.520503985230E*01 *.52507I?I6605E‘OI -.530854861470E“OI -.525663235152E‘OI -.498755067463E‘01 -.443733745582E’01 “.360955859047E-OI “.251251792681E-OI -.226889534304E-OI DIFFERENCE -.2185I4741391E~03 -.930888245435E‘O4 *.4095842541655‘04 -.248598487819E*03 -.2309I6855799E-03 -.I994I9I46160E‘04 -.10657OI43544E-03 v.453429037407E-03 -.584721099600E-03 -.337579377636E-03 -.254371625495E-04 -.10651I088235E-04 -.368635195635E-03 -.892378422577E-03 -.13361251013OE~02 ~.152453482682E-02 -.I42475150373E-02 -.104884464556E-02 -.95123513029BE-03 -.460899764228E~04 -.I97612425836E“03 -.727576789373E*04 -.489266803932E-04 -.246254222478E-03 -.164053064513E-03 .156308942678E*05 -.l94888150872E-03 -.492591637347E-03 -.449437306088E-03 -.139029818460E*03 .314305634682E-04 -.188268861169E-03 -.647269092049E-O3 -.IIII31679480E~02 -.133485571590E-02 -.128769335610E-02 -.951110498601E-03 -.858272491314E-03 K8110 Y/D .310 .820 .830 .340 .850 .860 .870 .880 .890 .900 .9I0 .920 .930 .940 .950 .960 .970 .980 .982 60 TABLE 6--CONTINUED PSI ‘.491299317750E“OI ‘.496833540572E-01 ”.50208933I203E‘OI “.506107013513E‘OI “.511997619338E‘OI ‘.516395901053E‘OI -.519187912083E*01 -.524676420639E‘OI -.529053012535E‘01 -.526936291033E‘01 ~.522959259382E“01 -.525058360916E‘OI ‘.53l48267386IE‘01 ~.530391057327E‘OI -.507855853444E-OI ".455636362247E-01 -.372563408215E-01 -.260020843823E*OI -.234770094637E-OI DIFFERENCE -.595893907317E-O4 -.531015202811E-04 -.179213l715816-03 -.337932669943E-04 -.8952912594545-04 -.23833562409oE-03 -.754313659869E-04 -.277o9345887oE-o4 -.3179I9594643E-03 -.468316521619E-03 -.245527414339E-03 .1335569322735-05 -.627812387407E-04 -.472782217003E-o3 -.910078598288E-03 -.ll9026166614E-02 -.116075491680E-02 -.87690511.1865-03 -.7aaoseoaaaaaa-oa II “ Illa-1|..IIII. l I: J.IIII ‘ o AT T=0 AFTER NOTE T TABLE 7 10 IBM 61 MOVEMENT 0F INTERFACE POSITION-‘CASE 5 (TABULATION FOR FIGURE 16) Y/D 0.800 0.810 0.820 0.830 0.840 0.850 0.860 0.870 0.880 0.890 0.900 0.910 0.920 0.930 0.940 0.950 0.960 0.970 0.980 0.982 .79570146E .80075022E .80438676E .80908718E .81980150E .82977453E .82755685E .83931947E .853363045 .82782176E .79248004E .81427781E .87429744E .85868187E .65180278E .21100515E -.42161320E -.11432038E “.91320570E .30086771E 1620 WAS USED 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 01 00 01 X/D .19211993E+02 .19290055E+02 .19364119E+02 .19434183E+02 .19500247E+02 .19562308E+02 .19620362E+02 .19674405E+02 .19724427E+02 .19770419E+02 .19812367E+02 .1985039IE+02 .19884412E+02 .19914431E+02 .19940447E+02 .19962461E+02 .19980472E+02 .19994481E+02 .20004487E+02 .20000000E+02 .19212544E 02 .19291271E 02 .19366286E 02 .19437256E 02 .19503400E 02 .19565676E 02 .19626150E 02 .19680301E 02 .19729976E 02 .1978456OE 02 .19839256E 02 .19876964E 02 .19898685E 02 .19939886E 02 .20060454E 02 .20331728E 02 .20843441E 02 .21754384E 02 .24120897E 02 .20903355E 02 FOR THE COMPUTATION. III I ‘ III-ll [III III nil II CHAPTER IV CONCLUSION This thesis has presented an analytical method by which the movement of the interface for transient two- phase flow in a porous medium can be predicted. Darcy's law is a well-established fundamental principle based upon which all problems concerning the flow through porous media are formulated. The use of Darcy's law in this con— tinuity equation then yields Laplace's equation, and the resulting boundary value problem is solved. Two approaches are used in solving the problem. The first of these is to consider the flow as one- dimensional while the other is to consider the flow as two-dimensional. The results of the onemdimensional analysis are represented by Figures 5 through 8. These results indi- cate that the velocity near the outflow seepage face cannot be approximated by the onemdimensional assumption, and the need for two—dimensional analysis is necessary. In the beginning of this investigation, the relaxation method was tried in solving the two-dimensional problem. However, it was found that this method would be unsuitable, if not impossible, for use with digital 62 63 computer, because the relaxation patterns in the vicinity of the interface change with the movement of the interface. This has led to the adoption of the present method. The results, as are represented by Figures 13 through 16, make available an analytical method which can be compared with field measurements and model tests. Although two dif- ferent sets of data are used in the two approaches, the results are qualitatively the same except near the out- flow seepage face. Comparison of Figures 8 and 16 reveals that while in the one-dimensional analysis reasonable results are obtained up to about-é : 2 from the end, a d reasonable result is obtained up to about-5 d = 0.1 in the two-dimensional analysis. The digital computer is utilized only for the two-dimensional analysis. The computation time for obtaining the results shown in Table 4 was about 60 minutes, while for Table 6 it was 160 minutes. For Tables 5 and 7, the time was just a fraction of a minute. It appears, therefore, that a computer which equals or excels Control Data 3600 in processing speed should be used. The method presented here can be applied to the problem with an initial interface of arbitrary shape. It can also be applied, with some modification, to the problem with different geometry. The problem in which the compressibility of fluids is taken into account and the effect of miscibility on the flow field present other possibilities for future study. l. 2. 3. 10. ll. BIBLIOGRAPHY Badon-Ghyben, W. Nota in verband met de voorgenomen putboring nabij Amsterdam, Ingenieur,gUtrecht, p. 21, 1888-89. Bear, J. and Dagon, G. Moving Interface in Coastal Aquifers, Journal ofiflydraulics Division, Vol. 90, No. HY4, July, 1964, p. 193. American Society of Civil Engineers. Collins, R. E. Flow of Fluids through Porous Materials, Reinhold, New York, 1961. Coats, K. H. Unsteady—State Liquid Flow through Porous Media Having Elliptic Boundaries, Petro- leum Transactions, AIME, T.N. 2050, Vol. 216, 1959. DeJong DeJosselin, G. Singularity Distributions for the Analysis of Multiple-Fluid Flow through Porous Media, Journal of Geo-physical Research, Vol. 65, No. 11, 1960. DeWiest, R. J. M. Unsteady-state Phenomena in Flow through Porous Media, Technical Report, No. 3, Department of Civil Engineering, Stanford Uni— versity, 1959. Dupuit, J. Etudes theoriques et pratiques sur le mouvementdes eaux a travers les terrains permeables, Carilian-Goeury, Paris, 1863. Forchheimer, Ph. Zeits Arch. Ing. Ver. Hannover, 1886. Henry, H. R. Salt Intrusion into Fresh-water Aquifers, Journal of Geophysical Research, Vol. 64, No. 11, 1959. Herzberg, B. Die Wasserversorgung einiger Nordse- ebader, J. Gasbeleuchtung and Wasserversorgung, 44, 815-19, 842-44, 1901. Hubbert, K. Darcy's Law and the Field Equations of the Flow of Underground Fluids, Journal of Petroleum Technology, October, 1956. 64 12. 13. 14. 15. 16. 17. 18. 19. 20. 65 Hubbert, M. K. The Theory of Ground-water Motion, Journal of Geology, 48, 785-944, 1940. Jacob, C. E. Flow of Ground Water, Chapter V, Engineering Hydraulics, edited by H. Rouse, John Wiley, New York, 1950. Kidder, R. E. Flow of Immiscible Fluids in Porous Media: Exact Solution of a Free Boundary Problem, Journal of Applied Physics, Vol. 27, No. 8, 1956. Kidder, R. E. Motion of the Interface between Two Immiscible Liquids of Unequal Density in a Porous Solid, Journal of Applied Physics, Vol. 27, No. 12, 1956. Kitagawa, K. Un aspect du develloppement des etudes des eaux soutteraines au Japon, Japanese Journal of Astronomy and Geophysics, 17, 141-55, 1939. Kunz, Kaiser S. Numerical Analysis, McGraw Hill, New York , 1957. Muskat, M. The Flow of Homogeneous Fluids through Porous Media, J. w. Edwards, Ann Arbor, 1946. Scheidegger, A. The Physics of Flow through Porous Media, Macmillan, New York, 1957. Todd, D. K. Sea Water Intrusion in Coastal Aquifers, Transactions, American Geophysical Union, Vol. 34, No. 5, 1953. APPENDIX I LIST OF SYMBOLS The following symbols are used throughout this 66 thesis: Symbol Definition Dimension al,a2 surface levels of liquid reservoir L b length of porous medium L d depth of porous medium L f(t) some function of time 9 acceleration due to gravity L/T2 h,h1,h2 piezometric head L 3 unit vector in the vertical direction K permeability of porous medium L2 k,kl,k2 hydraulic conductivity L/T kl product of hydraulic conductivity and buoyancy L/T m,n integers n1,n2 normals to a surface p pressure in fluid F/L2 Q discharge per unit width in porous 2 medium L /T If local velocity vector L/T 5 distance along the interface L t time T 67 Symbol Definition Dimension u,u1,u2 x components of‘E’ L/T v,vl,v2 y components of a. L/T x horizontal co~ordinate L Y vertical co-ordinate of interface L y vertical co-ordinate L '7, 71, 7% specific weights of fluids F/L3 F = 72 "' 2’1 3’1 /9 density of fluid M/L3 fl( 3 ,7 ) intensity of source f horizontal co-ordinate along interface L 7 vertical co~ordinate along interface L I? modified stream function L2/T 9V stream function L2/T ¢,¢1,¢2 velocity potentials L2/T 9 porosity of porous medium A! viscosity of fluid M/LT G( ) Green's function ( ) etc. partial derivatives of a function ( ) with respect to x, etc. Laplacian operator APPENDIX II FORMULAS USED IN NUMERICAL APPROXIMATION A. Newton's Divided-difference Interpolation Formula x u f(y) = f(y0) + (y - y0)f(y0,y1) + (y - y0)(y - yl) f(yo,y1,y2) + (y - yo)(y - y1)(y - y2)f(y0,yl,y2,y3) + . . . + (y‘- yo)(y - yl) . . . (y — yn_1)f(y0,y1, . . . yn) where f(yo,yl), f(y0,yl,y2) . . . etc. are the divided differences of f(y) defined as f(y0) - f(y1) f(Y $Y ) = _ O l y0 y1 f( ) f(y0.y1) - f(y1.y2) YO,Y1,Y2 = I _ y0 y2 f(yo.y1.y2. . . . yn> = f\Uo Then we have nnx mn um n(x,y) = cos ——— sin-~—x 1 b d and 2 2 A= -12 wry-5% b d 82 The function f(x,y) can be expanded in a double Fourier series, “3 d> n . f(x,y) = 2:1 220 (cos-BS25 Sin-agham’n where b d a =5;fofof(f,7)sinm—n-¥cosn—T%Edfd7. m,n The solution u can then be stated' thus: ‘“ a u(x,y) = - Z. Z gin 2 sing-Ex cos 2%. m=1 n=0 n 2025 +-E§) b d (On substituting in the above equation the expression for a , we obtain m,n 00 00 Z Z 4 cos __ngx sin Md b d f u(X,y) = - -77—- f( , ) ' m=O n=011bd (23+ELE) O O 7 b d nnrg . mniz cos b 5111 d <1de [b d 4 0° 9° - [—TX Z ° 0 O n bd m=0 n=0 cos ngx sin-9%x cos-32%: sin«EEgL 2 2 ]f(f,7)d}d7. n +‘E2 E7 d If the function of four variables in the square bracket be considered separately, ’Ibid., p. 69. by means of this function the solution of the problem can be very simply expressed: b d u(x,y) = f0 f0 G(x,y;},7)f(§,7)d}’d7. The function G(x,y: f, 7) is the Green's function for a rectangle.‘ B. Application to the Problem of Chapter III In the problem of Chapter III, the boundary portion L1 consists of the top and bottom of the rectangle of Figure 11, and the equation is Laplace's equation instead of Poisson's equation, so that the right-hand side f(x,y) is zero everywhere except at points on the interface. Since the interface is a line of sources, we have in effect taken it as the limiting case of a narrow zone, say of uniform width w on which f( E ,7) = f( {#7) where ’0( f, 7) is the source strength. Then the area integral ‘See pp. 384-86 in R. Courant and D. Hilbert, Methods of Mathematical Ph sics, Vol. I, Interscience PuSlIshers, New York, I953. 84 firea of G(x,y; f ,7 )f( f , ’Z)dA - Zone becomes the line integral . (3,7) jInterface [G(x,y,:§,>Z) w ~ ](wds) which in the limit (as w approaches zero) reduces to the line integral G(x,y; f, ’2) f( I, 7)ds to which must be added the integral along BA of Figure 11, the portion of L2 where the normal derivative is pres— cribed to be different from zero. The remaining part of and L will be zero in this case 1 2 because of the boundary conditions on 9?. the integrals along L Another way of looking at the problem is to note that each term of the series of Equation (50) satisfies Laplace's equation. Hence the series should satisfy Laplace's equation at interior points not on the inter— face, if it converges sufficiently. Each term of the series also satisfies the homogeneous boundary conditions. It remains only to show that the whole series satisfies the inhomogeneous boundary condition lpx = -k' on BA and the prescribed discontinuity condition on the interface. Since the Green's function is the sum of 3% ln-% and a function h(x,y) harmonic in the whole rectangle, the required normal derivative discontinuity can only be furnished by the logarithmic term. What we have to show in order to justify Equation (47) is that 85 B A _ .1 1 _ f .1 .1. f(x,y) — _[6 2n 1n r lfl(,§,’2)ds k B (Zn In r)dy has the prescribed normal derivative discontinuity, a jump of (0( f ,7) on DB and a jump of -k' across BA. This can be shown from the potential theory‘ which states that the potential of a simple distribution on a curve f(x,y) = .]C 7’ln-l ds r if the curve has a continuously turning tangent and 1f is bounded and integrable, is continuous for all finite points of the plane including passage through C. If C has continuous curvature and ‘7 is continuous, then the normal derivatives of f approach limits when P approaches A on C from either the positive or negative side, which satisfy the equations where A is point on the curve C. From this relation, f f ‘9 2 -.ii_l = -2n0£[. = - g 9 n1 9 n1 2% which is the required jump condition for f(x,y) along DB or BA. 'Quoted from W. J. Sternberg and T. L. Smith, The Theory of Potential and Spherical Harmonics, University of Toronto Press, 1944.