- Email: [email protected]
The finite element method: Is weighted volume integration essential?
Advances in Water Resources, Vol. 18, No. 6, pp. 345-352, 1995 Copyright 0 1995 Elsevier Science Lhited Printed in Great Britain. All rights reserved 0309-I 708/95 $09.50 + 0.00
0309-1708(95)00020-8
ELSEVIER
The finite element method: Is weighted volume integration essential? T. N. Narasimhan Materials Science and Mineral Engineering, 467 Evans Hall, University of California at Berkeley, Berkeley, California 94720-I 760, USA
(Received 5 December 1993; accepted 8 June 1995) ‘But the emperor has nothing on at all!!‘, said a little child. 1805-1875.
The Emperor’s New Clothes, Hans Christian Andersen,
In developing finite element equations for steady state and transient diffusiontype processes, weighted volume integration is generally assumed to be an intrinsic requirement. It is shown that such finite element equations can be developed directly and with ease on the basis of the elementary notion of a surface integral. Although weighted volume integration is mathematically correct, the algebraic equations stemming from it are no more informative than those derived directly on the basis of a surface integral. An interesting upshot is that the derivation based on surface integration does not require knowledge of a partial ditTerentia1 equation but yet is logically rigorous. It is commonly stated that weighted volume integration of the differential equation helps one carry out analyses of errors, convergence and existence, and therefore, weighted volume integration is preferable. It is suggested that because the direct derivation is logically consistent, numerical solutions emanating from it must be testable for accuracy and internal consistency in ways that the style of which may differ from the classical procedures of error- and convergence-analysis. In addition to simplifying the teaching of the finite element method, the thoughts presented in this paper may lead to establishing the finite element method independently in its own right, rather than it being a surrogate of the differential equation. The purpose of this paper is not to espouse any one particular way of formulating the finite element equations. Rather, it is one of introspection. The desire is to critically examine our traditional way of doing things and inquire whether alte:mate approaches may reveal to us new and interesting insights.
INTRODUCTION
equations. It appears that one can in fact successfully do so. An implication is that although the more elaborate derivations based on the Gale&in method or the variational principle are mathematically correct, they are no more informative or rigorous than the derivation based on a surface integral. The purpose of this paper is not to favor one particular way of deriving the finite element equations over another. Instead, the desire is to inquire whether by committing ourselves to the traditional way, we overlook useful alternate methods of solving problems. We begin with the derivation of the finite element equations for diffusion-type problems on the basis of surface integration, without invoking the partial differential equation. Following this we briefly examine how the finite element method has evolved over nearly four decades and reflect on how it may benefit from alternate perceptions.
It is generally assumed in the literature that weighted volume integration, either explicit (the Gale&n method; the method of weighted residuals) or implicit (minimization of the variational principle), is intrinsic to the development of the finite element method. In developing the finite element equations the common practice is to start with a weighted volume integral of the partial differential equation, apply Green’s identity to split the volume integral into a volume integral (for the interior domain) and a surface integral (for the boundary conditions) and use the resulting expression as the basis for writing down the discretized equations. The question we address in this communication is whether the very same finite element equations can be derived in a straight forward manner by starting with a simple surface integral and directly writing down the algebraic 345
T. N. Narasimhan
346
DERIVATION OF FINITE ELEMENT EQUATIONS The steady-state problem A basic premise of the steady-state diffusion process is that accumulation (mass, energy or other extensive quantity as appropriate) within any closed subdomain within the flow domain is zero. Using a surface integral, this can be expressed as
where q is flux density, n is the unit outer normal to the surface segment dI’ and I is the closed surface of the volume element V. Figure l(A) shows such an arbitrary subdomain. For computational convenience, let us consider a special case of Fig. l(A) in which the subdomain is a polygonal assemblage of triangles around a pointj, as in Fig. l(B). Let us use an equation of motion of the form, q = -KV$, where K is conductivity and 4 is the dependent variable which is a potential. Then, in view of eqn (l), the accumulation within the domain shown in Fig. l(B) is given by
where Aj is the accumulation in the domain around the pointj. The subscript (j, m) denotes the surface segment (that is, the side of the triangle) opposite nodal point, j, in triangular element m. For a steady-state system, this accumulation must be zero, according to eqn (1). Therefore, one could conceivably start with a linear interpolation function for 4 defined over a triangle, derive the algebraic relations relevant to the discretized surface integral (2) using elementary principles of analytical geometry, set it to zero and solve for the potentials. How will this procedure compare with the conventional finite element procedure? Recall that the accumulation term is evaluated in the conventional finite element approach with the help of the weighted volume integral, A;=21
Vcj m=l
* KV&$k
dV
vrn
where AT is the accumulation around j obtained by the conventional finite element procedure involving weighted volume integration, cj is the weighting function with respect to nodal point j and the subscript k denotes the three nodal points of element m. It has been shown by previous workers (e.g. Wilson;” Narasimhan and Witherspoon”) that the algebraic reduction of the weighted volume integral (3) leads to
CJ 5
V
m=l Vrn =;
(fl ux entering the outer boundary of domain)
= i 2 KV$j,, *nj,,AI’j,, m=l
-
A.
Note that the last summation in eqn (4) is precisely the same as the discretized surface integral (2) of the surface integral (l), except for the multiplier l/2. Thus, we see that the weighted volume integral yields an accumulation Ai which is exactly one half of the surface integral (Aj) evaluated over the entire domain (Fig. l(B)) surrounding j. That is, Aj = *(AT). If we have a steady-state problem, setting Aj = 0 or setting Ai = 0 and solving the algebraic equations should lead to identical results. The transient problem
Fig. 1. Domain around nodal point j defined by its neighbors: (A) an arbitrarily shaped element, and (B) special case of an assemblage of triangular elements.
Let us now consider the transient problem. In this case the accumulation must be absorbed within the domain of integration, leading to a change in the dependent variable 4 at the nodal points as dictated by the capacitance of the materials within the domain. Consider, for convenience, one of the triangular subdomains of Fig. l(B), as shown in Fig. 2. The flux entering this element across the side opposite j must be partitioned
Is weighted volume integration essential?
347
form (Felippa’),
VCj*KVkh d V
kl A.
B.
Fig. 2. (A) Equal partitioning of the capacitance of a triangular
element among its comer points, obtained by drawing medians, and (B) partitioning of flux among the three corner points. and absorbed by the three nodal points by virtue of their capacitances. Intuitively one would suspect that the partitioning will depend on the nature of the local flow field over the time interval of interest. As the flow field changes in time, so will the relative amounts partitioned between the nodes. However, if we do not assume prior knowledge of such a flow field, we need to predetermine a reasonable logic for partitioning the accumulation. Accordingly, if we show no bias to any of the three nodes, then each of them will have to share the element equally in regard to capacitance, as shown by the dotted lines in Fig. 2(A), obtained by joining the medians. Now consider the tot,al flux F entering the surface (j, m) as shown in Fig. 2(B). As suggested in the figure, exactly one half of this total flux enters the domain controlled byj and the remaining half is shared without bias by the other two nodes. The half which enters the domain of j can be expressed as being equal to (1/3)Vec’A$j/At, where A4j is the change in potential at nodej. Here, Ve is the volume of the element and ce is the specific capacitance of the material contained in the element. Now since the remaining half is shared equally by the other two nodes, these accumulations will, respectively, be equal to, (l/6) V’c”A&i /At and (l/6) VeceAq&/At. Consequently, we obtain the mass balance equation for node j on the basis of surface integration,
Aj = - e
=
qj,m * nj,,AI’j,,
1 (5)
where k,,, and k2,m are the other two nodes of element m besides j. Now let us comparce this with the discretized equation of the conventional finite element method in the context of the distributed capacitance matrix
If we recognize that Aj G 2Ai, it readily follows that eqns (5) and (6) are fully equivalent. That is, the derivation based on the simple surface integral leads exactly to the same discretized equation as that obtained by weighted volume integration or by minimizing the variational principle. Note that eqn (5) does not constitute the only logical way to partition the accumulation. In arriving at eqn (5) we saw that half of the accumulation is accommodated by node j and the remaining half is accommodated collectively by the nodes other than j. Therefore, one may equally well write two logically valid mass balance statements as follows. A.
A=__
2
1’ c 2,,1
qj,m *nj,mAQm =
and,
a,_ --- l5c qj,m *nj,mArj,m 2
2m=1
=
1
The expression in eqn (7) is used by many in the literature instead of eqn (6) as the basis for the finite element method. This expression is sometimes referred to as the ‘lumped’ capacitance matrix (Neuman and Narasimhani2). Note particularly that eqn (8) gives rise to a capacitance matrix which involves only off-diagonal terms. Logically, eqn (8) appears to be as credible as eqns (7) or (5) and hence could form the basis for solving for the dependent variable. Yet, to the knowledge of this writer, this form of the capacitance matrix has not been recognized so far in the literature. Whereas the weights of the off-diagonal terms are lumped into the diagonal term to obtain eqn (7) the weight of the diagonal term is distributed into the off-diagonal term in eqn (8). One could reasonably speculate that eqns (7) and (8) should lead to closely similar, if not identical solutions. The right hand side of eqn (8) may be called a ‘deleted-diagonal’ capacitance matrix.
DISCUSSION
If indeed this simple derivation from a surface integral yields the same result as that stemming from the more
T. N. Narasimhan
abstract mathematical derivations, why is it that this line of thought has not received attention in the past from previous workers. Are there other advantages inherent in the conventional derivation of the finite element method which are lost when we go directly from the surface integral to the algebraic equations? Or indeed, is it possible that the simpler derivation offers hitherto unrecognized advantages? We proceed to answer these questions by first briefly recapitulating the evolution of ideas relating to the finite element method. Evolution of ideas As is well known, the finite element method originated in the structural engineering literature. Dating back to the 192Os, ideas pertaining to the method have historically developed along two parallel lines. Practicing structural engineers visualized large complex structures as assemblages of numerous structural elements of simple shapes and the stiffness of the complete structure was obtained by summing the stiffness of the individual components. In essence, based on physics and intuition, the engineering approach involved direct integration of discrete quantities. The turn of the century also saw the application of variational calculus to the solution of problems of mathematical physics through the seminal contributions of Rayleigh and Ritz. Consequently, mathematicians such as Courant who were exposed to the engineering solution of structural problems, saw the methods of Rayleigh and Ritz latent in the physically based methodology of the structural engineers. In retrospect, we find each of these approaches to be self consistent and mutually equivalent. Yet, when we look at modern literature on the finite element method dealing with diffusion-type problems, we find that the mathematical approach is given far greater attention than the intuitive engineering approach. The phrase ‘finite element method’ was first used by Clough.3 In this classic paper, Clough3 derived stiffness of typical structural elements of triangular and rectangular shapes and showed how the stiffnesses of these discrete components can be summed up subject to conditions of force equilibrium and compatibility of deformation. He further demonstrated that the resulting set of linear algebraic equations can be advantageously solved with the help of the digital computer. It is worth noting here that Clough derived the equations on the basis of physical reasoning rather than by mathematically integrating the differential equation. The finite element concept as embodied in Clough3 evolved out of an earlier paper in the aeronautics literature by Turner et aZ.,14 which too relied on physical reasoning in developing the stiffness matrix. Indeed, these authors referred to their technique as the ‘direct stiffness method’. If we recognize that the term stiffness pertains to a discrete structural component of prescribed shape, size and
material make-up, we readily see that these early workers directly wrote down the discrete integral expressions from physical principles rather than going through an intermediate phase of integrating a differential equation. Although the finite element method crystallized in its present form in the late 1950s rudimentary notions relevant to the method seem to have existed much earlier. According to Finlayson and Scriven,6 Biezeno’ used a particular form of weighted volume integration, the subdomain method, to solve for stress-strain equilibrium in two-dimensional media. In doing so, Biezeno was drawing upon the earlier ideas of Ritz and Galerkin. Biezeno’s paper contains a comment by Courant pointing out that the subdomain method could be viewed as a logical consequence of variational calculus. Subsequently, in an address delivered before the American Mathematical Society in 1941, Courant4 elaborated on how variational methods can be applied to solve problems or equilibrium and vibrations. Some researchers (R.L. Taylor, personal communication, 1994) would consider that Courant’s4 paper constitutes the first formal effort at providing a mathematical foundation for the finite element method. Wilson” was among the early workers to apply the finite element method to transient diffusion-type problems. In the process, he too used physical reasoning to set up the ‘lumped’ capacitance matrix as shown in eqn (7). Unlike the equilibrium problems of stress analysis for which the earliest finite element method was devised, the nature of the diffusion process is such that a true variational principle does not exist for the transient problem. To overcome this, Gurtin’ proposed the use of a variational principle defined in the Laplace domain. Subsequent workers (e.g. Javandel and Witherspoon’) took advantage of Gurtin’s work in developing finite element equations by the RayleighRitz method. In parallel, the Galerkin method (e.g. Pinder and Frind13) was also being used for developing the equations. Thus, as we look at the early developments leading to the finite element method, we find that sound physical reasoning as well as abstract mathematical ideas have played parallel roles. A desire for mathematical rigor motivated many researchers to choose the variational approach and the Galerkin procedure in preference to the intuitive approach. One of the early consequences of this preference was the introduction of isoparametric elements. However, not everyone was completely persuaded by the purely mathematical development of the equations. Numerical difficulties encountered in solving non-linear problems and multi-phase problems often necessitated the moderation of mathematical formalism with physical intuition. For example, Neuman” encountered numerical difficulties in solving the non-linear unsaturated groundwater flow problem when using the
Is weighted volume integration essential?
distributed capacitance matrix. To overcome this, he used physical reasoning and resorted successfully to the use of a lumped capacitance matrix. Even as rigorous mathematics was being used to firmly establish the finite element method, some workers (e.g. Narasimhan and Witherspoon;” Neuman and Narasimhan12) pointed out the advantages of following a physically-based approach. As the finite element method was being applied increasingly to non-linear diffusion-type problems in the earth sciences and in petroleum engineering, more and lmore researchers have resorted to combine intuition with mathematical formalism. This is evidenced by the fact that the 1980s has seen a growing popularity of what has been called the controlvolume finite element method, which is essentially a variant of the approach pioneered by Wilson” nearly two decades earlier. With this brief background on how ideas have evolved with reference t’o the finite element method, let us now revert to a discussion of surface integration in comparison with weigh&d volume integration. Surface integration
and weighted volume integration
We begin by restricting consideration to the robust linear basis function (triangular element in two dimensions). The equiva,lence between eqns (5) and (6) suggests that the information content of the two discretized equations are exactly the same. Thus it is clear that no special advantage is to be gained in the weighted volume integration procedure, except for the fact that the capacitance weighting was automatically taken care of by the mlathematical process without a need for physical judgement: we have, as it were, a default option built into weighted volume integration to generate capacitance welights. One may speculate as to why derivation based on the surface integral is not preferred in the finite element literature. Perhaps the answer lies in a tradition of mathematical physics. It is that dynamic systems are to be described in terms of one or more differential equations. Although the power of integral formulations has been recognized in science from around the turn of the century (Cajori2), the tendency to think almost exclusively in terms of the differential equations is so strong that the paradigm of numerical modeling is that they are approximate solvers of differential equations. Therefore, it is not surprising that many researchers prefer to base the finite element method on the partial differential equation. Nevertheless, the fact that we can logically derive the finite element equations independent of the differential equation suggests that the current paradigm is open to question. Perhaps the direct derivation deserves a little more attention. Although functional analysis provides for error-, convergence- and! existence-analysis, it does not preclude the possibility that other approaches could be
349
&P A.
B.
Fig. 3. Triangular finite elements in two different types of local
patterns of flow geometry, (A) uniform flow, (B) convergent flow. developed to evaluate the credibility of solutions generated on the basis of surface integration. As an example, let us examine how we may go about evaluating the accuracy of the surface integration procedure leading to the accumulation term A, in eqn (2). Consider a triangular element shown in Fig. 3(A). As shown, this subdomain exists in a region where the flow field is locally uniform and parallel, at the scale of discretization. Clearly, in this case, linear interpolation over the triangular domain is an excellent approximation, errors will be minimal and accuracy high. In Fig. 3(B), the element exists in a region where the flow field is converging and non-uniform on the scale of discretization. In this case linear interpolation over the triangle is clearly not an accurate assumption. How well may we quantitatively estimate the accuracy of a solution based on the linear basis function? Not very well, for the following reason. In general, convergent and divergent flow fields will arise in diffusion-type problems due to a variety of causes such as complex boundary shapes and heterogeneity. In transient problems, the flow patterns may even change as a function of time. Hence, magnitude of local numerical errors will be highly problem-dependent as well as time-dependent. Therefore, one cannot carry out universal error analysis or convergence analysis valid for all problems. This line of reasoning is provided merely to suggest that functional analysis does not constitute the only means of estimating accuracy. Other styles of evaluating internal consistency and accuracy of numerical solutions should be possible. In order that errors may be minimized, one will have to refine the scale of discretization in those areas where non-uniform flow patterns exist, based on the premise that at sufficiently fine scales all flows can be reasonably approximated as uniform and parallel. Extension
to isoparametric
elements
Conceptually, the higher-order isoparametric elements are no different from the simpler triangular elements for the steady-state case. Figure 4 shows a node j at the intersection of four isoparametric elements. In this case too, surface integration over the closed surface r will yield the accumulation within the domain. Under steady
T. N. Narasimhan
350
Q
A
A
”
A.
B.
Fig. 5. Two ways of partitioning an isoparametric element, (A)
physically unrealistic, and (B) physically more realistic.
Fig. 4. Domain around nodal point, j, defined by four
isoparametric elements. state, by setting this integral to zero and solving the algebraic equation must yield the same results that would be obtained if one used a conventional weighted volume integration procedure. As a matter of detail, it may be noticed here that if one chooses to follow surface integration, then one can dispense with Gaussian quadrature within the element and replace it with quadrature at appropriate points on I to accurately evaluate the surface fluxes. Indeed, Gaussian quadrature on the surface may entail fewer integration points than is the case with conventional Galerkin procedure, thus contributing to computational efficiency while still preserving accuracy. Suppose the problem involving the higher order elements is transient. In this case, use of the surface integration procedure would require that an appropriate logic is provided to partition the accumulation among the various nodal points. Unfortunately, this is a difficult task. Unlike the triangular element, in which the partitioning with the help of medians was physically credible, an isoparametric element with nodes only on the sides of the element, is difficult to partition in a physically meaningful way among the nodal points defining the element. For example, if one takes a central point and partitions the domain as shown in Fig. 5(A), the partitioning is clearly unrealistic. Instead, if one partitions the region reasonably as shown in Fig. 5(B), the shaded region in the middle is unaccounted for in the mass balance evaluation. One way out may be to use the partition in Fig. 5(B) and distribute the accumulation in the shaded region among the other nodes, weighted according to the volume of the partition commanded by each node.
Thus, in the case of the isoparametric element under transient conditions, one may consider the physically based approach to be inadequate. At the same time the weighted volume integration scheme is particularly attractive because the capacitance weights are automatically provided by weighting functions. Caution, however, is in order. The Galerkin weights are purely intended for spatial interpolation. No compelling rationale exists to assert that the spatial weight is indeed appropriate for weighting capacitances in the transient problem. As pointed out earlier, the actual weights of integration will be closely related to local flow geometry and hence will be a function of both space and time. Therefore, a set of capacitance weights chosen a priori to be independent of space and time will not necessarily be adequate, although it may superficially circumvent a troublesome issue. Under the circumstances, it may in fact be preferable to look for weights guided by a physical understanding of the particular problem than to rely on a default option. Thus, it may well be that our inability to partition the isoparametric element in a satisfactory way is an indication that the isoparametric element is logically less than adequate for solving transient problems. This may be the reason why many researchers still use triangular elements in preference to isoparametric elements for transient problems. Intrinsic advantage of the finite element concept At this juncture it is worth asking, ‘What attributes of the finite element concept set it clearly apart from the classical finite differences? The answer is that the finite element method enables the simultaneous evaluation of spatial gradients in more than one arbitrarily chosen direction at a given location in the flow domain, whereas classical finite differences could evaluate gradients only in orthogonal directions. It is not surprising that the finite element notion originated in the stress analysis literature. To simultaneously evaluate normal stress and shear stresses or, equivalently, normal strains and distortion at a given location, one has to use tensorial quantities. In turn, tensorial quantities have to be used
Is weighted volume integration
in conjunction with potential gradients normal to and parallel to surface segments of interest. When the domain of interest has non-trivial geometry, the surface segments may not everywhere be oriented parallel to a simple Cartesian system. Thus, structural engineers dealing with complex shapes such as air-foils imaginatively devised the finite element logic to provide them the necessary flexibility of evaluating spatial gradients in any arbitrary direction. It is important to recognize here that the procedure of weighted volume integration is not essential to take advantage of this special ability of the finite element method. The power of the finite element logic can be exploited quite well with the surface integration approach. Accuracy and credibility of solutions In formulating the finite element equations, we are concerned with issues such as rigor, accuracy and model credibility. A few thoughts are presented below to suggest how one may assert the credibility of a solution when the finite element equations are derived with the help of the surface integral, without resorting to a differential equation. We start by recognizing that in diffusion-type problems fluxes are related to spatial variations of potentials through empirically established equations of motion. In general, fluxes are inversely related to resistances. Equation (1) implies that in a steady-state system the algebraic sum of fluxes over any closed subdomain must be zero. Since flux is related to resistance it follows that if the resistances between nodes are calculated correctly in the physically derived finite element method, then the solution generated will automatically be accurate. Suppose we are using linear basis functions over triangular finite elements. If the mesh is coarse and the actual variation of potential deviates from linearity, then the error associated with the deviation is incorporated into the stiffness matrix, which, in diffusion-type problems is indeed a resistance matrix. Consequently, one can argue that instead of following the usual procedure of error analysis, one could equally well evaluate the accuracy of finite element solutions by focussing attention on how accurately resistances are being calculated. Noting that resistances are governed by geometry of flow patterns in addition to material properties, we recognize that a need exists to develop new approaches to evaluate the credibility of numerical solutions by focussing attention on flow geometry. There is reason to believe that this approach, unlike the conventional error analysis approach, will have the ability to evaluate the credibility of any solution generated by the finite element model, however complex the problem may be. For transient problems, accuracy of calculated resistances has to be supplemented by the accuracy of calculated capacitances. If, in the final matrices which
essential?
351
are used in the calculation of changes in potential with time, one can assert that resistances and capacitances are correct, then the solution generated with these matrices must be correct; one does not have to seek a proof based on functional analysis.
CONCLUDING
REMARKS
As can be seen from the literature, the development of the finite element equations for flow and transport in groundwater systems is almost entirely based (explicitly or implicitly) on weighted volume integration. Nevertheless, it appears that the same equations can be developed in a straightforward manner by starting from a surface integral and directly writing down the algebraic equations. An interesting upshot of such a development is that the notion of a differential equation does not enter the picture at all. If indeed the reasoning presented in this paper is correct, the implication is that we will have a finite element formulation which is rigorous but yet is independent of the differential equation and the variational principle. Such an independent base should provide the motivation to search for alternate ways of addressing issues related to model verification and error analysis which are different in style from the conventional wisdom based on the theory of differential equations. Over 25 years ago Finlayson and Scriven7 pointed out that a variational principle is really not necessary if one is interested merely in formulating the discretized finite element equations; one could achieve the same end much more easily through weighted volume integration. The reasoning presented in this paper suggests that even weighted volume integration is not really needed to formulate the finite element equations; a surface integral is more than adequate for the purpose.
ACKNOWLEDGEMENTS Thanks are due to W. G. Gray, Pierre Perrochet and Andrew Tompson for a review of the manuscript.
REFERENCES Biezeno, C. B., Graphical and numerical methods for solving stress problems. In Proc. 1st Int. Cong. for Appl. Mech., De& 1924. Delft-Technische Boekhandel en Drukkerij J. Waltman Jr, 1925, pp. 3-17. Cajori, F., A History of Mathematics (3rd edn). Chelsea Publishing Company, New York, USA, 1980. Clough, R. W., The finite element method in plane stress analysis. In 2nd Conf. Electronic Computation. Am. Sot. Civil Engng, 1960, pp. 345-77. Courant, R., Variational methods for the solution of problems of equilibrium and vibrations. Bull. Am. Math. Sot., 49 (1943) l-23.
T. N. Narasimhan 5. Felippa, C. A., Refined Finite Element Analysis for Linear and Non-Linear Two Dimensional Structures. Report No. SESM 66-22, Structural Engineering Laboratory, Department of Civil Engineering, University of California at Berkeley, CA, USA, 1966, 6. Finlayson, B. A. and Striven, L. E., The method of weighted residuals-a review. Appl. Me&. Rev., 19(9) (1966) 735-48. 7. Finlayson, B. A. and Striven, L. E., On the search for variational principles. J. Heat Muss Trans., 10 (1967) 799. 8. Gurtin, M. E., Variational principles for linear initialvalue problems. Quart. Appl. Math., 22(3) (1964) 252-6. 9. Javandel, I. and Witherspoon, P. A., Application of the finite element method to transient flow in porous media. Sot. Petrol Engng J., 8(3) (1968) 241-52. 10. Narasimhan, T. N. and Witherspoon, P. A., An integrated finite difference method for fluid flow in porous media. Water Resow. Res., 12(l) (1976) 57-64.
11. Neuman, S. P., Saturated-unsaturated seepage by finite elements. J. Hydruul. Div., Am. Sot. Civil Engng, 99 (1975) 2233-50. 12. Neuman, S. P. and Narasimhan, T. N., Mixed explicitimplicit finite element scheme for diffusion-type problems Part I: Theory. Znt. J. Num. Meths Engng, ll(2) (1977) 30924. 13. Pinder, G. F. and Frind, E. O., Application of Galerkin procedure to aquifer analysis. Water Resow. Res., 8 (1972) 108-20. 14. Turner, M. J., Clough, R. W., Martin, H. C. and Topp, L. J., Stiffness and deflection analysis of complex structures. J. Aeronaut. Sci., 23(9) (1956) 805-23. 15. Wilson, E. L., The Determination of Temperatures Within Mass Concrete Structures. Report No. SESM 68-17, Structural Engineering Laboratory, Department of Civil Engineering, University of California at Berkeley, CA, USA, 1968.
0309-1708(95)00020-8
ELSEVIER
The finite element method: Is weighted volume integration essential? T. N. Narasimhan Materials Science and Mineral Engineering, 467 Evans Hall, University of California at Berkeley, Berkeley, California 94720-I 760, USA
(Received 5 December 1993; accepted 8 June 1995) ‘But the emperor has nothing on at all!!‘, said a little child. 1805-1875.
The Emperor’s New Clothes, Hans Christian Andersen,
In developing finite element equations for steady state and transient diffusiontype processes, weighted volume integration is generally assumed to be an intrinsic requirement. It is shown that such finite element equations can be developed directly and with ease on the basis of the elementary notion of a surface integral. Although weighted volume integration is mathematically correct, the algebraic equations stemming from it are no more informative than those derived directly on the basis of a surface integral. An interesting upshot is that the derivation based on surface integration does not require knowledge of a partial ditTerentia1 equation but yet is logically rigorous. It is commonly stated that weighted volume integration of the differential equation helps one carry out analyses of errors, convergence and existence, and therefore, weighted volume integration is preferable. It is suggested that because the direct derivation is logically consistent, numerical solutions emanating from it must be testable for accuracy and internal consistency in ways that the style of which may differ from the classical procedures of error- and convergence-analysis. In addition to simplifying the teaching of the finite element method, the thoughts presented in this paper may lead to establishing the finite element method independently in its own right, rather than it being a surrogate of the differential equation. The purpose of this paper is not to espouse any one particular way of formulating the finite element equations. Rather, it is one of introspection. The desire is to critically examine our traditional way of doing things and inquire whether alte:mate approaches may reveal to us new and interesting insights.
INTRODUCTION
equations. It appears that one can in fact successfully do so. An implication is that although the more elaborate derivations based on the Gale&in method or the variational principle are mathematically correct, they are no more informative or rigorous than the derivation based on a surface integral. The purpose of this paper is not to favor one particular way of deriving the finite element equations over another. Instead, the desire is to inquire whether by committing ourselves to the traditional way, we overlook useful alternate methods of solving problems. We begin with the derivation of the finite element equations for diffusion-type problems on the basis of surface integration, without invoking the partial differential equation. Following this we briefly examine how the finite element method has evolved over nearly four decades and reflect on how it may benefit from alternate perceptions.
It is generally assumed in the literature that weighted volume integration, either explicit (the Gale&n method; the method of weighted residuals) or implicit (minimization of the variational principle), is intrinsic to the development of the finite element method. In developing the finite element equations the common practice is to start with a weighted volume integral of the partial differential equation, apply Green’s identity to split the volume integral into a volume integral (for the interior domain) and a surface integral (for the boundary conditions) and use the resulting expression as the basis for writing down the discretized equations. The question we address in this communication is whether the very same finite element equations can be derived in a straight forward manner by starting with a simple surface integral and directly writing down the algebraic 345
T. N. Narasimhan
346
DERIVATION OF FINITE ELEMENT EQUATIONS The steady-state problem A basic premise of the steady-state diffusion process is that accumulation (mass, energy or other extensive quantity as appropriate) within any closed subdomain within the flow domain is zero. Using a surface integral, this can be expressed as
where q is flux density, n is the unit outer normal to the surface segment dI’ and I is the closed surface of the volume element V. Figure l(A) shows such an arbitrary subdomain. For computational convenience, let us consider a special case of Fig. l(A) in which the subdomain is a polygonal assemblage of triangles around a pointj, as in Fig. l(B). Let us use an equation of motion of the form, q = -KV$, where K is conductivity and 4 is the dependent variable which is a potential. Then, in view of eqn (l), the accumulation within the domain shown in Fig. l(B) is given by
where Aj is the accumulation in the domain around the pointj. The subscript (j, m) denotes the surface segment (that is, the side of the triangle) opposite nodal point, j, in triangular element m. For a steady-state system, this accumulation must be zero, according to eqn (1). Therefore, one could conceivably start with a linear interpolation function for 4 defined over a triangle, derive the algebraic relations relevant to the discretized surface integral (2) using elementary principles of analytical geometry, set it to zero and solve for the potentials. How will this procedure compare with the conventional finite element procedure? Recall that the accumulation term is evaluated in the conventional finite element approach with the help of the weighted volume integral, A;=21
Vcj m=l
* KV&$k
dV
vrn
where AT is the accumulation around j obtained by the conventional finite element procedure involving weighted volume integration, cj is the weighting function with respect to nodal point j and the subscript k denotes the three nodal points of element m. It has been shown by previous workers (e.g. Wilson;” Narasimhan and Witherspoon”) that the algebraic reduction of the weighted volume integral (3) leads to
CJ 5
V
m=l Vrn =;
(fl ux entering the outer boundary of domain)
= i 2 KV$j,, *nj,,AI’j,, m=l
-
A.
Note that the last summation in eqn (4) is precisely the same as the discretized surface integral (2) of the surface integral (l), except for the multiplier l/2. Thus, we see that the weighted volume integral yields an accumulation Ai which is exactly one half of the surface integral (Aj) evaluated over the entire domain (Fig. l(B)) surrounding j. That is, Aj = *(AT). If we have a steady-state problem, setting Aj = 0 or setting Ai = 0 and solving the algebraic equations should lead to identical results. The transient problem
Fig. 1. Domain around nodal point j defined by its neighbors: (A) an arbitrarily shaped element, and (B) special case of an assemblage of triangular elements.
Let us now consider the transient problem. In this case the accumulation must be absorbed within the domain of integration, leading to a change in the dependent variable 4 at the nodal points as dictated by the capacitance of the materials within the domain. Consider, for convenience, one of the triangular subdomains of Fig. l(B), as shown in Fig. 2. The flux entering this element across the side opposite j must be partitioned
Is weighted volume integration essential?
347
form (Felippa’),
VCj*KVkh d V
kl A.
B.
Fig. 2. (A) Equal partitioning of the capacitance of a triangular
element among its comer points, obtained by drawing medians, and (B) partitioning of flux among the three corner points. and absorbed by the three nodal points by virtue of their capacitances. Intuitively one would suspect that the partitioning will depend on the nature of the local flow field over the time interval of interest. As the flow field changes in time, so will the relative amounts partitioned between the nodes. However, if we do not assume prior knowledge of such a flow field, we need to predetermine a reasonable logic for partitioning the accumulation. Accordingly, if we show no bias to any of the three nodes, then each of them will have to share the element equally in regard to capacitance, as shown by the dotted lines in Fig. 2(A), obtained by joining the medians. Now consider the tot,al flux F entering the surface (j, m) as shown in Fig. 2(B). As suggested in the figure, exactly one half of this total flux enters the domain controlled byj and the remaining half is shared without bias by the other two nodes. The half which enters the domain of j can be expressed as being equal to (1/3)Vec’A$j/At, where A4j is the change in potential at nodej. Here, Ve is the volume of the element and ce is the specific capacitance of the material contained in the element. Now since the remaining half is shared equally by the other two nodes, these accumulations will, respectively, be equal to, (l/6) V’c”A&i /At and (l/6) VeceAq&/At. Consequently, we obtain the mass balance equation for node j on the basis of surface integration,
Aj = - e
=
qj,m * nj,,AI’j,,
1 (5)
where k,,, and k2,m are the other two nodes of element m besides j. Now let us comparce this with the discretized equation of the conventional finite element method in the context of the distributed capacitance matrix
If we recognize that Aj G 2Ai, it readily follows that eqns (5) and (6) are fully equivalent. That is, the derivation based on the simple surface integral leads exactly to the same discretized equation as that obtained by weighted volume integration or by minimizing the variational principle. Note that eqn (5) does not constitute the only logical way to partition the accumulation. In arriving at eqn (5) we saw that half of the accumulation is accommodated by node j and the remaining half is accommodated collectively by the nodes other than j. Therefore, one may equally well write two logically valid mass balance statements as follows. A.
A=__
2
1’ c 2,,1
qj,m *nj,mAQm =
and,
a,_ --- l5c qj,m *nj,mArj,m 2
2m=1
=
1
The expression in eqn (7) is used by many in the literature instead of eqn (6) as the basis for the finite element method. This expression is sometimes referred to as the ‘lumped’ capacitance matrix (Neuman and Narasimhani2). Note particularly that eqn (8) gives rise to a capacitance matrix which involves only off-diagonal terms. Logically, eqn (8) appears to be as credible as eqns (7) or (5) and hence could form the basis for solving for the dependent variable. Yet, to the knowledge of this writer, this form of the capacitance matrix has not been recognized so far in the literature. Whereas the weights of the off-diagonal terms are lumped into the diagonal term to obtain eqn (7) the weight of the diagonal term is distributed into the off-diagonal term in eqn (8). One could reasonably speculate that eqns (7) and (8) should lead to closely similar, if not identical solutions. The right hand side of eqn (8) may be called a ‘deleted-diagonal’ capacitance matrix.
DISCUSSION
If indeed this simple derivation from a surface integral yields the same result as that stemming from the more
T. N. Narasimhan
abstract mathematical derivations, why is it that this line of thought has not received attention in the past from previous workers. Are there other advantages inherent in the conventional derivation of the finite element method which are lost when we go directly from the surface integral to the algebraic equations? Or indeed, is it possible that the simpler derivation offers hitherto unrecognized advantages? We proceed to answer these questions by first briefly recapitulating the evolution of ideas relating to the finite element method. Evolution of ideas As is well known, the finite element method originated in the structural engineering literature. Dating back to the 192Os, ideas pertaining to the method have historically developed along two parallel lines. Practicing structural engineers visualized large complex structures as assemblages of numerous structural elements of simple shapes and the stiffness of the complete structure was obtained by summing the stiffness of the individual components. In essence, based on physics and intuition, the engineering approach involved direct integration of discrete quantities. The turn of the century also saw the application of variational calculus to the solution of problems of mathematical physics through the seminal contributions of Rayleigh and Ritz. Consequently, mathematicians such as Courant who were exposed to the engineering solution of structural problems, saw the methods of Rayleigh and Ritz latent in the physically based methodology of the structural engineers. In retrospect, we find each of these approaches to be self consistent and mutually equivalent. Yet, when we look at modern literature on the finite element method dealing with diffusion-type problems, we find that the mathematical approach is given far greater attention than the intuitive engineering approach. The phrase ‘finite element method’ was first used by Clough.3 In this classic paper, Clough3 derived stiffness of typical structural elements of triangular and rectangular shapes and showed how the stiffnesses of these discrete components can be summed up subject to conditions of force equilibrium and compatibility of deformation. He further demonstrated that the resulting set of linear algebraic equations can be advantageously solved with the help of the digital computer. It is worth noting here that Clough derived the equations on the basis of physical reasoning rather than by mathematically integrating the differential equation. The finite element concept as embodied in Clough3 evolved out of an earlier paper in the aeronautics literature by Turner et aZ.,14 which too relied on physical reasoning in developing the stiffness matrix. Indeed, these authors referred to their technique as the ‘direct stiffness method’. If we recognize that the term stiffness pertains to a discrete structural component of prescribed shape, size and
material make-up, we readily see that these early workers directly wrote down the discrete integral expressions from physical principles rather than going through an intermediate phase of integrating a differential equation. Although the finite element method crystallized in its present form in the late 1950s rudimentary notions relevant to the method seem to have existed much earlier. According to Finlayson and Scriven,6 Biezeno’ used a particular form of weighted volume integration, the subdomain method, to solve for stress-strain equilibrium in two-dimensional media. In doing so, Biezeno was drawing upon the earlier ideas of Ritz and Galerkin. Biezeno’s paper contains a comment by Courant pointing out that the subdomain method could be viewed as a logical consequence of variational calculus. Subsequently, in an address delivered before the American Mathematical Society in 1941, Courant4 elaborated on how variational methods can be applied to solve problems or equilibrium and vibrations. Some researchers (R.L. Taylor, personal communication, 1994) would consider that Courant’s4 paper constitutes the first formal effort at providing a mathematical foundation for the finite element method. Wilson” was among the early workers to apply the finite element method to transient diffusion-type problems. In the process, he too used physical reasoning to set up the ‘lumped’ capacitance matrix as shown in eqn (7). Unlike the equilibrium problems of stress analysis for which the earliest finite element method was devised, the nature of the diffusion process is such that a true variational principle does not exist for the transient problem. To overcome this, Gurtin’ proposed the use of a variational principle defined in the Laplace domain. Subsequent workers (e.g. Javandel and Witherspoon’) took advantage of Gurtin’s work in developing finite element equations by the RayleighRitz method. In parallel, the Galerkin method (e.g. Pinder and Frind13) was also being used for developing the equations. Thus, as we look at the early developments leading to the finite element method, we find that sound physical reasoning as well as abstract mathematical ideas have played parallel roles. A desire for mathematical rigor motivated many researchers to choose the variational approach and the Galerkin procedure in preference to the intuitive approach. One of the early consequences of this preference was the introduction of isoparametric elements. However, not everyone was completely persuaded by the purely mathematical development of the equations. Numerical difficulties encountered in solving non-linear problems and multi-phase problems often necessitated the moderation of mathematical formalism with physical intuition. For example, Neuman” encountered numerical difficulties in solving the non-linear unsaturated groundwater flow problem when using the
Is weighted volume integration essential?
distributed capacitance matrix. To overcome this, he used physical reasoning and resorted successfully to the use of a lumped capacitance matrix. Even as rigorous mathematics was being used to firmly establish the finite element method, some workers (e.g. Narasimhan and Witherspoon;” Neuman and Narasimhan12) pointed out the advantages of following a physically-based approach. As the finite element method was being applied increasingly to non-linear diffusion-type problems in the earth sciences and in petroleum engineering, more and lmore researchers have resorted to combine intuition with mathematical formalism. This is evidenced by the fact that the 1980s has seen a growing popularity of what has been called the controlvolume finite element method, which is essentially a variant of the approach pioneered by Wilson” nearly two decades earlier. With this brief background on how ideas have evolved with reference t’o the finite element method, let us now revert to a discussion of surface integration in comparison with weigh&d volume integration. Surface integration
and weighted volume integration
We begin by restricting consideration to the robust linear basis function (triangular element in two dimensions). The equiva,lence between eqns (5) and (6) suggests that the information content of the two discretized equations are exactly the same. Thus it is clear that no special advantage is to be gained in the weighted volume integration procedure, except for the fact that the capacitance weighting was automatically taken care of by the mlathematical process without a need for physical judgement: we have, as it were, a default option built into weighted volume integration to generate capacitance welights. One may speculate as to why derivation based on the surface integral is not preferred in the finite element literature. Perhaps the answer lies in a tradition of mathematical physics. It is that dynamic systems are to be described in terms of one or more differential equations. Although the power of integral formulations has been recognized in science from around the turn of the century (Cajori2), the tendency to think almost exclusively in terms of the differential equations is so strong that the paradigm of numerical modeling is that they are approximate solvers of differential equations. Therefore, it is not surprising that many researchers prefer to base the finite element method on the partial differential equation. Nevertheless, the fact that we can logically derive the finite element equations independent of the differential equation suggests that the current paradigm is open to question. Perhaps the direct derivation deserves a little more attention. Although functional analysis provides for error-, convergence- and! existence-analysis, it does not preclude the possibility that other approaches could be
349
&P A.
B.
Fig. 3. Triangular finite elements in two different types of local
patterns of flow geometry, (A) uniform flow, (B) convergent flow. developed to evaluate the credibility of solutions generated on the basis of surface integration. As an example, let us examine how we may go about evaluating the accuracy of the surface integration procedure leading to the accumulation term A, in eqn (2). Consider a triangular element shown in Fig. 3(A). As shown, this subdomain exists in a region where the flow field is locally uniform and parallel, at the scale of discretization. Clearly, in this case, linear interpolation over the triangular domain is an excellent approximation, errors will be minimal and accuracy high. In Fig. 3(B), the element exists in a region where the flow field is converging and non-uniform on the scale of discretization. In this case linear interpolation over the triangle is clearly not an accurate assumption. How well may we quantitatively estimate the accuracy of a solution based on the linear basis function? Not very well, for the following reason. In general, convergent and divergent flow fields will arise in diffusion-type problems due to a variety of causes such as complex boundary shapes and heterogeneity. In transient problems, the flow patterns may even change as a function of time. Hence, magnitude of local numerical errors will be highly problem-dependent as well as time-dependent. Therefore, one cannot carry out universal error analysis or convergence analysis valid for all problems. This line of reasoning is provided merely to suggest that functional analysis does not constitute the only means of estimating accuracy. Other styles of evaluating internal consistency and accuracy of numerical solutions should be possible. In order that errors may be minimized, one will have to refine the scale of discretization in those areas where non-uniform flow patterns exist, based on the premise that at sufficiently fine scales all flows can be reasonably approximated as uniform and parallel. Extension
to isoparametric
elements
Conceptually, the higher-order isoparametric elements are no different from the simpler triangular elements for the steady-state case. Figure 4 shows a node j at the intersection of four isoparametric elements. In this case too, surface integration over the closed surface r will yield the accumulation within the domain. Under steady
T. N. Narasimhan
350
Q
A
A
”
A.
B.
Fig. 5. Two ways of partitioning an isoparametric element, (A)
physically unrealistic, and (B) physically more realistic.
Fig. 4. Domain around nodal point, j, defined by four
isoparametric elements. state, by setting this integral to zero and solving the algebraic equation must yield the same results that would be obtained if one used a conventional weighted volume integration procedure. As a matter of detail, it may be noticed here that if one chooses to follow surface integration, then one can dispense with Gaussian quadrature within the element and replace it with quadrature at appropriate points on I to accurately evaluate the surface fluxes. Indeed, Gaussian quadrature on the surface may entail fewer integration points than is the case with conventional Galerkin procedure, thus contributing to computational efficiency while still preserving accuracy. Suppose the problem involving the higher order elements is transient. In this case, use of the surface integration procedure would require that an appropriate logic is provided to partition the accumulation among the various nodal points. Unfortunately, this is a difficult task. Unlike the triangular element, in which the partitioning with the help of medians was physically credible, an isoparametric element with nodes only on the sides of the element, is difficult to partition in a physically meaningful way among the nodal points defining the element. For example, if one takes a central point and partitions the domain as shown in Fig. 5(A), the partitioning is clearly unrealistic. Instead, if one partitions the region reasonably as shown in Fig. 5(B), the shaded region in the middle is unaccounted for in the mass balance evaluation. One way out may be to use the partition in Fig. 5(B) and distribute the accumulation in the shaded region among the other nodes, weighted according to the volume of the partition commanded by each node.
Thus, in the case of the isoparametric element under transient conditions, one may consider the physically based approach to be inadequate. At the same time the weighted volume integration scheme is particularly attractive because the capacitance weights are automatically provided by weighting functions. Caution, however, is in order. The Galerkin weights are purely intended for spatial interpolation. No compelling rationale exists to assert that the spatial weight is indeed appropriate for weighting capacitances in the transient problem. As pointed out earlier, the actual weights of integration will be closely related to local flow geometry and hence will be a function of both space and time. Therefore, a set of capacitance weights chosen a priori to be independent of space and time will not necessarily be adequate, although it may superficially circumvent a troublesome issue. Under the circumstances, it may in fact be preferable to look for weights guided by a physical understanding of the particular problem than to rely on a default option. Thus, it may well be that our inability to partition the isoparametric element in a satisfactory way is an indication that the isoparametric element is logically less than adequate for solving transient problems. This may be the reason why many researchers still use triangular elements in preference to isoparametric elements for transient problems. Intrinsic advantage of the finite element concept At this juncture it is worth asking, ‘What attributes of the finite element concept set it clearly apart from the classical finite differences? The answer is that the finite element method enables the simultaneous evaluation of spatial gradients in more than one arbitrarily chosen direction at a given location in the flow domain, whereas classical finite differences could evaluate gradients only in orthogonal directions. It is not surprising that the finite element notion originated in the stress analysis literature. To simultaneously evaluate normal stress and shear stresses or, equivalently, normal strains and distortion at a given location, one has to use tensorial quantities. In turn, tensorial quantities have to be used
Is weighted volume integration
in conjunction with potential gradients normal to and parallel to surface segments of interest. When the domain of interest has non-trivial geometry, the surface segments may not everywhere be oriented parallel to a simple Cartesian system. Thus, structural engineers dealing with complex shapes such as air-foils imaginatively devised the finite element logic to provide them the necessary flexibility of evaluating spatial gradients in any arbitrary direction. It is important to recognize here that the procedure of weighted volume integration is not essential to take advantage of this special ability of the finite element method. The power of the finite element logic can be exploited quite well with the surface integration approach. Accuracy and credibility of solutions In formulating the finite element equations, we are concerned with issues such as rigor, accuracy and model credibility. A few thoughts are presented below to suggest how one may assert the credibility of a solution when the finite element equations are derived with the help of the surface integral, without resorting to a differential equation. We start by recognizing that in diffusion-type problems fluxes are related to spatial variations of potentials through empirically established equations of motion. In general, fluxes are inversely related to resistances. Equation (1) implies that in a steady-state system the algebraic sum of fluxes over any closed subdomain must be zero. Since flux is related to resistance it follows that if the resistances between nodes are calculated correctly in the physically derived finite element method, then the solution generated will automatically be accurate. Suppose we are using linear basis functions over triangular finite elements. If the mesh is coarse and the actual variation of potential deviates from linearity, then the error associated with the deviation is incorporated into the stiffness matrix, which, in diffusion-type problems is indeed a resistance matrix. Consequently, one can argue that instead of following the usual procedure of error analysis, one could equally well evaluate the accuracy of finite element solutions by focussing attention on how accurately resistances are being calculated. Noting that resistances are governed by geometry of flow patterns in addition to material properties, we recognize that a need exists to develop new approaches to evaluate the credibility of numerical solutions by focussing attention on flow geometry. There is reason to believe that this approach, unlike the conventional error analysis approach, will have the ability to evaluate the credibility of any solution generated by the finite element model, however complex the problem may be. For transient problems, accuracy of calculated resistances has to be supplemented by the accuracy of calculated capacitances. If, in the final matrices which
essential?
351
are used in the calculation of changes in potential with time, one can assert that resistances and capacitances are correct, then the solution generated with these matrices must be correct; one does not have to seek a proof based on functional analysis.
CONCLUDING
REMARKS
As can be seen from the literature, the development of the finite element equations for flow and transport in groundwater systems is almost entirely based (explicitly or implicitly) on weighted volume integration. Nevertheless, it appears that the same equations can be developed in a straightforward manner by starting from a surface integral and directly writing down the algebraic equations. An interesting upshot of such a development is that the notion of a differential equation does not enter the picture at all. If indeed the reasoning presented in this paper is correct, the implication is that we will have a finite element formulation which is rigorous but yet is independent of the differential equation and the variational principle. Such an independent base should provide the motivation to search for alternate ways of addressing issues related to model verification and error analysis which are different in style from the conventional wisdom based on the theory of differential equations. Over 25 years ago Finlayson and Scriven7 pointed out that a variational principle is really not necessary if one is interested merely in formulating the discretized finite element equations; one could achieve the same end much more easily through weighted volume integration. The reasoning presented in this paper suggests that even weighted volume integration is not really needed to formulate the finite element equations; a surface integral is more than adequate for the purpose.
ACKNOWLEDGEMENTS Thanks are due to W. G. Gray, Pierre Perrochet and Andrew Tompson for a review of the manuscript.
REFERENCES Biezeno, C. B., Graphical and numerical methods for solving stress problems. In Proc. 1st Int. Cong. for Appl. Mech., De& 1924. Delft-Technische Boekhandel en Drukkerij J. Waltman Jr, 1925, pp. 3-17. Cajori, F., A History of Mathematics (3rd edn). Chelsea Publishing Company, New York, USA, 1980. Clough, R. W., The finite element method in plane stress analysis. In 2nd Conf. Electronic Computation. Am. Sot. Civil Engng, 1960, pp. 345-77. Courant, R., Variational methods for the solution of problems of equilibrium and vibrations. Bull. Am. Math. Sot., 49 (1943) l-23.
T. N. Narasimhan 5. Felippa, C. A., Refined Finite Element Analysis for Linear and Non-Linear Two Dimensional Structures. Report No. SESM 66-22, Structural Engineering Laboratory, Department of Civil Engineering, University of California at Berkeley, CA, USA, 1966, 6. Finlayson, B. A. and Striven, L. E., The method of weighted residuals-a review. Appl. Me&. Rev., 19(9) (1966) 735-48. 7. Finlayson, B. A. and Striven, L. E., On the search for variational principles. J. Heat Muss Trans., 10 (1967) 799. 8. Gurtin, M. E., Variational principles for linear initialvalue problems. Quart. Appl. Math., 22(3) (1964) 252-6. 9. Javandel, I. and Witherspoon, P. A., Application of the finite element method to transient flow in porous media. Sot. Petrol Engng J., 8(3) (1968) 241-52. 10. Narasimhan, T. N. and Witherspoon, P. A., An integrated finite difference method for fluid flow in porous media. Water Resow. Res., 12(l) (1976) 57-64.
11. Neuman, S. P., Saturated-unsaturated seepage by finite elements. J. Hydruul. Div., Am. Sot. Civil Engng, 99 (1975) 2233-50. 12. Neuman, S. P. and Narasimhan, T. N., Mixed explicitimplicit finite element scheme for diffusion-type problems Part I: Theory. Znt. J. Num. Meths Engng, ll(2) (1977) 30924. 13. Pinder, G. F. and Frind, E. O., Application of Galerkin procedure to aquifer analysis. Water Resow. Res., 8 (1972) 108-20. 14. Turner, M. J., Clough, R. W., Martin, H. C. and Topp, L. J., Stiffness and deflection analysis of complex structures. J. Aeronaut. Sci., 23(9) (1956) 805-23. 15. Wilson, E. L., The Determination of Temperatures Within Mass Concrete Structures. Report No. SESM 68-17, Structural Engineering Laboratory, Department of Civil Engineering, University of California at Berkeley, CA, USA, 1968.