Zienkiewicz's Work

During his years at Edinburgh, he largely devoted himself to work on lubrication theory. An important outcome of this work is described in Reference 2, in which a significant new theory of lubrication is advanced, demonstrating that lift can be obtained due to a temperature distribution across the film. Moving to Northwestern, he devoted himself principally to topics in solid mechanics. A striking, and indeed landmark paper, was Reference 3 which combines linear programming and plastic design theory into a unified computer-oriented approach to the elastic design of beam and frame structures. A. Charnes, his co-author, is a major figure in the development of linear programming and the paper describes a general approach to the problem and sets out, perhaps for the first time, the connection between duality in linear programming algorithms and duality in the formulation of the limit analysis problem. It is one of the earliest papers to utilise the new methodology of mathematical programming in structural mechanics.

Another series of papers from the Northwestern era deserving of particular attention are those written with R. W. Gerstner. The set of papers is mainly concerned with a general approach to elasticity problems (e.g. Reference 4) in which exact and numerical solutions are combined, where the numerical solution is done as a correction to the exact solution.

It should be recognised that when Olek and his colleagues began research into the finite element method at Swansea in 1961, the method was already a known technique. The mathematical foundations of the technique were laid in the early decades of the last century, with contributions by Ritz, Galerkin, Courant and others. Within the finite element community, the origins of the method are attributed to the contribution of Hrenikoff & McHenry (1941) who, using engineering intuition, approximated the response of a plate by an equivalent framework system and Courant (1943) who developed a more rigorous mathematical approach to the solution of a torsion problem by subdivision of the domain into triangular sub-regions. Progress from that time onward was slow, no doubt hampered by the lack of computing power, until the first recognisable finite element paper, by Turner, Clough, Martin & Topp appeared in 1956. By the end of that decade a substantial body of literature had been accumulated and the term ‘finite element’ was coined by R. W. Clough, University of California at Berkeley in 1960. It was against this background that Olek entered the finite element research arena. It should also be mentioned that it was at this time that J. H. Argyris FRS, Imperial College and later the University of Stuttgart, embarked on finite element research and their careers ran in parallel for the next four decades, both contributing prominently to the field. To use Zienkiewicz’s own words (5) to describe this scene:

“My own introduction to finite elements dates back to 1958 when I first met Clough. With much of my earlier work being based on finite difference calculus (and relaxation methods – as an ex-member of the Southwell team in the 1940s), I naturally contested at first the merits of the new method which seemed to concentrate on the solution of problems which by then we were fully capable of solving. However, in one area of structural analysis the merit of subdivision seemed obvious. That was the field of shell analysis where the difficulties of formulating differential equations in curvilinear coordinates and subsequently applying finite difference approximations could be easily sidestepped by considering a shell to be composed of a series of simple, flat, triangular elements. It was in this area of activity that I made my first start (1962) – not yet realising that the search for a good triangular element to represent plate flexure was to be fraught with difficulties. This problem also motivated the research of Ray Clough and with a friendly competition we arrived at our goal some three years later presenting a series of triangular plate and shell elements at the Wright-Patterson Conference of 1965”

The Wright-Patterson Conference held in Dayton, Ohio in 1965 was a landmark in finite element analysis and paper (6) describing the first triangular plate bending elements from the Swansea group contained a remarkable number of major new ideas. One was the use of area coordinates for integration over triangular regions. Another of the contributions in the paper was the notion of the patch test. The paper had described two finite elements, one of which was non-conforming. In a displacement finite element formulation, displacement continuity is assured across element boundaries, due to description of the element displacement field in terms of (common) interfacial nodal displacements, but the interface element tractions are generally discontinuous since the element stresses are derived from the, approximated, displacement field on an individual element basis. Elements of this type are termed as being non-conforming. Such elements were not new in finite element analysis, but no one had attempted to establish criteria for their convergence.

Bruce Irons, then working at Rolls Royce and co-author of the paper – and who was destined to become a colleague at Swansea, identified the conditions to be met in the choice of element approximation functions and, in the case that the condition of the inter-element conformity was not met, devised a simple test to be applied to a collection, or ‘patch’, of elements. The patch test has proven to be of fundamental importance in the finite element theory. It was also at this meeting that J. Tinsley Oden, a pioneer of the finite element method, made first contact with Olek (7):

“I first met Olek at the Dayton finite element meetings in the mid 1960s, which is where some say the finite element method, being born in the mid 1950s, reached its adolescence. There were a number of original and important works that formed the foundations of the subject that were presented there by engineers and scientists working in this new and exciting field. Of course, Olek was already known to many there because of his first textbook on finite elements, co-authored with Y. K. Cheung. Olek’s intense interest and warmth was intriguing. We hit it off immediately and began a friendship that lasted until his passing some four decades later”

Work on the notion of the patch test continued over the years. In collaboration with Robert L. Taylor, he extended the procedure to mixed element formulations, in which both displacement and stress terms are considered as primary variables (8) and which was subsequently used to ensure convergence of some new element forms for plates (9, 10).

It was during this period that the first industrial application of the finite element method, in Europe at least, took place in 1963 when Zienkiewicz and co-workers undertook the stress analysis of the Clywedog dam in mid-Wales. As can be seen from Figure 2, which illustrates the dam and the finite element discretisation employed, the mesh was extremely coarse by present day standards; which reflected the limited computational power available at that time. Nevertheless, the analysis provided valuable insight to the behaviour of the dam and its foundation.

A paper that appeared in 1965 and which was to have profound impact in later years was ‘Finite Elements in the Solution of Field Problems’ (11), co-authored with Y. K. Cheung. Up to this time the method had been restricted to structural mechanics problems, by expanding techniques for the analysis of frameworks into a method for the analysis of elastic continua. As such, the methodology relied heavily of the theorem of total potential energy. Zienkiewicz was able to perceive it as a more general tool for the analysis of all types of problems in mathematical physics that could be described in terms of a differential equation with given boundary and initial conditions. By employing weighted residual procedures, and in particular a Galerkin approach, computational solutions could be obtained for classes of problem where a potential function cannot be easily derived. In particular, he identified the procedure as a scheme for solving Laplace’s equation which governs the behaviour of ideal fluids and the torsion of prismatic sections. Olek and his colleagues soon amplified these ideas to deal with problems of heat transfer (12) and electromagnetics (13).

It is well established that isoparametric representation is fundamental to many modern aspects of finite element analysis. The essence of the idea originated in contributions by Ian Taig of the British Aircraft Corporation (a forerunner of British Aerospace) and Bruce Irons, then at Rolls Royce. Taig described a scheme for distorting the shape of a rectangular element to quadrilateral form in a contribution to the 1963 Swansea Symposium (14). Subsequently, Irons published a note in the AIAA Journal which described, in generality, the approach as it is known today and coined the term ‘isoparametric’.

The contribution of Reference 15 from Swansea brought together the ideas put forward in the earlier contributions and formalised the procedure for general application. It was around this time that Thomas J. R. Hughes first saw Olek in person at a seminar at the Massachusetts Institute of Technology (16).

“There was a lot of excitement in the room because Olek was about to present the ‘isoparametric concept’, which ultimately revolutionised the subject, but the image that remains in my mind was he began the lecture with the picture of a man made out of finite elements. It got a big laugh. At the time it was so incongruous to think of a man made out of finite elements, but now it is a serious business. Who could have thought at that time we would be eventually using medical imaging data to build anatomically and physiologically realistic models of real people for simulating medical interventions”.

This early picture, intended to convey the message that the study of finite elements should be an enjoyable pursuit, is reproduced in Figure 3.

The use of flat plate finite elements for the analysis of curved shell structures, although satisfactory for a number of situations, has severe limitations in others. This became clear towards the end of the 1960s, motivating a widespread activity in the formulation of curved shell finite elements. The obvious avenue towards this goal was via the established thin shell theory. Zienkiewicz, however, had other ideas and had long expressed his interest in the analysis of shells by a simpler approach. In place of invoking the Love-Kirchhoff theory, he advocated dispensing with any shell theory and using general finite elements. The idea was rendered especially attractive on account of the emergence of the concept of isoparametric representation, which would facilitate the description of the curved geometry.

The first elaboration of the notion of degenerating solids to create shell elements appeared in1968 (17) at the Second Conference on Matrix Methods in Structural Mechanics at Wright-Patterson Air Force Base, Ohio and somewhat later in 1970 (18). These papers set the stage for exploration of avenues which were subsequently pursued by many authors. However, the behaviour of truly thin shells could not be replicated. The breakthrough in overcoming this situation came in 1971 when Zienkiewicz, R. L. Taylor and J. M. Too (19) identified the source of difficulty as being an unrealistic representation of transverse shear phenomena. The difficulty was resolved through the use of a reduced numerical integration quadrature in evaluating the shear contribution to the stiffness term which was of a lower order than would be required for exact evaluation of this integral. This operation was termed ‘reduced integration’. Bob Taylor describes the background to this work as follows (20):

“In 1968 when travelling to the Second Conference on Matrix Methods in Structural Mechanics in Dayton, Ohio, I met Olek for the first time in person and I was pleased to be invited to Swansea for the academic year 1969-70. Fate stepped in that year to keep Olek near Swansea as he had an operation to adjust the alignment of his leg to ease difficulties in walking. During his post-operative recovery in Singleton Hospital and later at his home, we met regularly to discuss the possibility of solving shell problems using curved three-dimensional elements – often over a glass of whiskey, which unfortunately may have contributed to a blood clot in his leg that prolonged his recovery period. Working with Olek’s doctoral student James Too, we did succeed in developing a successful shell formulation based on an isoparametric brick combined with shell approximations and reduced integration”.

The reduced integration concept found application, in later years, in work in other aspects of the finite element method. A somewhat different direction in plate bending analysis was taken in 1978 (21) where the Reissner-Mindlin theory was adopted, which involved independent description of the transverse displacement and the slope (rotation) fields. Again, the reduced integration concept is essential to the capturing of thin plate behaviour. By this time, however, the order of the integration to be employed could be related on a rational basis to the geometry of the element and the assumed approximating functions.

In 1966 Olek and his colleagues turned their attention to inelastic analysis. Time dependant problems were the first in the class of inelastic problems to draw Zienkiewicz’s attention. In 1966, with M. Watson, he formulated finite element procedures for the analysis of concrete structures undergoing creep deformation (22). Then, with Watson and King (23), he considered the general viscoelastic problem within the framework of finite element analysis.

Reference (24) represented a major contribution to the finite element analysis of time-independent inelastic problems. To be sure, contributions by others had preceded this paper, but this work brought together the characterisation of material behaviour in terms of incremental theory and a new procedure for advancing along the load-displacement curve. Termed the ‘initial stress’ approach it represented, at that time, a fundamental advancement in the computational modelling of inelastic problems.

Collaborations with two Ph.D. students resulted in further developments in the field of inelastic analysis. Extensive studies in the modelling of plasticity problems were undertaken with G. C. Nayak (25). Three noteworthy contributions are the so-called ‘alpha-constant stiffness’ algorithm for iterative solution of the algebraic equations arising in nonlinear finite element analysis (26), a convenient form of the third invariant of the stress tensor suitable for yield surface computations (27) and a review which unified the constitutive relations for elasto-plastic finite element analysis (28).

The other collaboration, with I. Cormeau, involved the concept of viscoplasticity. Olek had for a long while envisioned a scheme whereby both time-dependent and time-independent inelastic phenomena could be handled by the same algorithmic structure. The first paper developing this approach appeared in 1972 (29). A by-product of this effort was the identification of a rational criterion for the determination of the time-step size for instability in time-dependent nonlinear analyses. The viscoplastic approach later found extensive usage in the analysis of extrusion and forming processes (30) and a range of geotechnical problems (31, 32).

Around this time the concept of an ‘overlay model’ was introduced for materials where characterisation of the constitutive response in terms of a single yield surface is inadequate. The overlay approach, in which a series of yield surfaces is introduced, was formulated by Zienkiewicz and co-workers in (33, 34).

The work in flow problems began in the 1965 paper with Cheung, as discussed earlier and, continued on to the treatment of seepage, electromagnetics and heat transfer. This set the stage for consideration of the more sophisticated problems which involve convective terms and nonlinearities, such as viscous, incompressible flow described by the Navier-Stokes equations. Zienkiewicz’s contributions to this latter topic began to appear in 1973, especially in 1974 in his Introductory Address to the First Symposium on Finite Elements in Flow Problems (35).

However, substantial problems remained for the solution of practical problems. One such difficulty, which was already well recognised in the finite difference literature, related to the fact that in the momentum equations, written in terms of velocities, the convection terms are proportional to the first derivative while the viscous term is proportional to the second derivative. The latter dominates at low values of the Reynolds number while the reverse is true at high Reynolds number. There is numerical instability as one increases from low to high speeds and to cope with this a technique known as ‘upwinding’ had been devised by finite difference researchers. 
Adaptation of this process to upwind weighting of the Petrov-Galerkin kind in the context of the finite element approximation was at first unclear. After much effort, Olek was able to identify the proper approach, and following numerical experimentation and theoretical refinement, a procedure was described in 1977 in the paper with Heinrich et al. (36) permitting the solution of both simple convective problems and of the incompressible Navier Stokes equations over a large range of Reynolds numbers. This development of upwinding in a finite element setting provided a general theoretical basis that did not previously exist.

However, this first step allowing the treatment of steady state problems was not sufficient, as the use of the upwind Petrov-Galerkin procedure was not justified when source or transient terms existed. Here, a major step forward was taken by the introduction of a characteristic Galerkin process making use of the wave propagation features of these together with the optimality of Galerkin methods for the self adjoint, diffusion part of the problem. This led to a new family of algorithms (37) which permitted a wide and new range of problems to be solved. The algorithm produced, in which the characteristic ‘search’ is dealt with by a simple Taylor expansion, could be reinterpreted as a finite element version of the Lax-Wendroff method and also identified with the so-called Taylor-Galerkin method. The latter can deal with a number of characteristic velocities and, as such, allowed extension to problems of high speed compressible gas flow to be made by the mid 1980s.

Another major difficulty in the numerical analysis of viscous incompressible flow problems arises on account of the constraint equations resulting from the incompressibility condition. Although there were conventional means of handling such constraint equations, e.g. through the use of Lagrange multipliers, a new and possibly more convenient approach to accommodating these constraints by means of a penalty method was used for the first time in 1974 in finite element fluid flow analysis by Godbole and Zienkiewicz (38). Successful applications, refinements and extensions of the penalty approach were subsequently reported by Zienkiewicz and his colleagues for the plastic flow of metals, assumed stress formulation of the finite element method, heat transfer, shallow water equations, natural convection and viscous flow analysis.

The group at Swansea contributed to the finite element analysis of dynamic problems from its earliest years in the field, with such developments as the consistent mass matrix (39) and the work of Irons on condensation schemes. A generalisation of the work in structural dynamics took place in 1973 (40), where Zienkiewicz demonstrated that certain of the existing time-stepping algorithms could be obtained by application of weighted residual methods to the equations of motion. In addition, some algorithms that had not previously been identified could be constructed in this way. This work is part of the more general notion of ‘finite elements in time’ (41). With W. L. Wood and R. L. Taylor (42) he showed that the weighting of independently interpolated displacement and velocity vectors, which is in correspondence with the well-known Newmark method under certain circumstances, allows the selection of time steps purely on the basis of structural response.

Geotechnical engineering was the focus of Olek’s attention for much of his professional career. His earlier work in this area was restricted to rock mass problems and his interest in soil mechanics only began in the 1970s. In reference (43) with D. J. Naylor, he adapted the popular, modern ‘critical state’ description of soil behaviour to finite element analysis. Except for cyclic load cases, this proved to be an effective model. Viscoplastic and more traditional inelastic constitutive relationships and finite element algorithms were applied in Reference 31 and, with G. N. Pande (44) he explored a number of candidate forms of yield surfaces for both soil and rock mechanics. Liquefaction continued to concern him and in Reference (45) a new theory was presented by him and co-workers. Many of the problems related to wave action and geotechnical foundations involve boundary conditions at infinity. Zienkiewicz and his colleagues, most notably P. Bettess, created and applied families of ‘infinite elements’ (46, 47) to account for such conditions.

Zienkiewicz also contributed to many other areas of computational modelling, but space only allows a brief mention of some of the more interesting topics. Early in his career at Swansea, he collaborated with B. Fraeijs de Veubeke in a treatment of plate bending formulations based on the complementary principle and assumed stress function fields (48). Zienkiewicz and R. L. Taylor returned to this theme in 1979 when they presented a method based on assumed stress fields and utilised the penalty function approach to enforce essential constraint conditions. During the 1970s Olek was intrigued by the connection between the boundary integral and finite element methods. He presented comprehensive studies (49, 50) examining the relationship, and was able to place them under a unified framework.