<figure><img src="image_1.png" /></figure>

HAL
open science

# A general methodology for recovering equilibrating finite element tractions and stress fields for plate and solid elements

Pierre Ladevèze, Edward Anthony Ward Maunder

► To cite this version:

Pierre Ladevèze, Edward Anthony Ward Maunder. A general methodology for recovering equilibrating finite element tractions and stress fields for plate and solid elements. Computer Assisted Methods in Engineering and Science, 1997, 4 (3-4), pp.533-548. hal-05223732

HAL Id: hal-05223732

https://hal.science/hal-05223732v1

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A general methodology for recovering equilibrating
finite element tractions and stress fields for plate
and solid elements

Pierre Ladevèze

Laboratoire de Mecanique et Technologie, E.N.S. Cachan / C.N.R.S./Universite Paris 6,
61 Avenue du President Wilson, 94235 Cachan Cedex, France

Edward Anthony Ward Maunder

School of Engineering, University of Exeter, Exeter, EX4 4QF, England

(Received December 20, 1996)

A recent geometric presentation of a general and efficient methodology for recovering equilibrating tractions and stress fields from 2-D conforming displacement finite element models is reviewed, and further considered in the context of plate elements in 2-D and 3-D. This methodology requires the resolution of corner nodal forces/moments, and this presents localised problems which are solved in a simple way by exploiting Maxwell force diagrams. For more complex 3-D models including solid elements where higher degrees of element connectivity occur, the geometric procedure is adapted so as to retain the computational efficiency gained from recognising the topological and geometrical properties of a finite element model. Graph-theoretic and algebraic topological concepts are invoked in this context. The equilibrating tractions recovered for each element enable statically admissible stresses to be computed element by element, and local Trefftz fields may be exploited.

**1. INTRODUCTION**

Conventional finite element models of structures are based on displacement fields which are gen-
erally conforming, but the corresponding stress fields are only in equilibrium in an “average” or
weak sense. Lack of equilibrium can lead to disastrous consequences, e.g. the destruction of the
Sleipner offshore platform has been partly attributed to serious underestimation of shear forces
based on finite element results [6]. Weak equilibrium may be strengthened as deemed necessary by
an appropriate adaptive procedure involving h- and/or p-refinement. The designer/analyst may re-
quire equilibrating solutions at various structural levels e.g. as resultant forces on sections, traction
distributions on sections, substructures, or individual elements; or as local stress fields. Various
methods have been proposed for recovering stronger forms of equilibrium from conventional models
[2, 9, 10, 11, 18] with varying degrees of complexity and success. The procedure [11] which recovers
equilibrium directly from a displacement model without further refinement was initially proposed
in the context of error estimation. Recent collaboration has led to a new interpretation and presen-
tation of this procedure, and an extension from regular to irregular meshes of membrane elements
[14]. In the present paper this new “geometrical” presentation is reviewed for the 2-D membrane
models in Section 2, before extending it in Section 3 to assemblies of plate elements where moment
vectors as well as force vectors exist in the modelling of plates, folded plates, or shells in 3-D.
Section 4 then extends the geometrical concepts of the methodology to plate and solid elements,
where higher degrees of element connectivity exist. Graph-theoretic concepts are here exploited in
new procedures aimed at reducing the computational cost of implementing the methodology.
## 2. REVIEW OF CONCEPTS FROM 2-D MEMBRANE ELEMENTS

Analysis of a conforming displacement model results in stresses $\sigma_h$ due to the prescribed loads. These stresses and the loads form the data for a postprocessing procedure to recover equilibrating tractions and stresses. The procedure can be summarised in the following five stages:

### (i) Identify nodal forces from the weak form of equilibrium of the finite element stress tensor $\sigma_h$

For a model with domain $\Omega$ and boundary $\Gamma$, $\sigma_h$ satisfies the integral Equation (1):

$$ - \int_{\Omega} (\sigma_h \cdot \operatorname{grad} \phi_n) \, d\Omega + \int_{\Omega} (\mathbf{f} \cdot \phi_n) \, d\Omega + \int_{\Gamma_t} (\mathbf{t} \cdot \phi_n) \, d\Gamma = 0, \quad (1) $$

where $\phi_n$ denotes the scalar shape function for node $n$, and $\mathbf{f}$ and $\mathbf{t}$ denote prescribed body forces on $\Omega$ and tractions on $\Gamma_t$ respectively. The third integral is present only when node $n$ lies on the boundary $\Gamma_t$. Each domain integral is the sum of contributions from element integrals which are interpreted as element nodal forces associated with node $n$ of element $E$:

$$ \mathbf{F}_n^E = \int_{\Omega^E} (\sigma_h^E \cdot \operatorname{grad} \phi_n^E) \, d\Omega; \quad \mathbf{f}_n^E = \int_{\Omega^E} (\mathbf{f}^E \cdot \phi_n^E) \, d\Omega; \quad \mathbf{T}_n^E = \int_{\Gamma_t^E} (\mathbf{t}^E \cdot \phi_n^E) \, d\Gamma, \quad (2) $$

$\mathbf{F}_n^E$ represents the response of an element to nodal displacements and is normally available in the output from a conventional analysis, $\mathbf{f}_n^E$, $\mathbf{T}_n^E$ are the consistent nodal loads which may be supplied as data to the displacement model. Equation (1) can now be interpreted as expressing the conditions for equilibrium of node $n$ as a free body subject to interactions with the elements $E$ connected to it and the consistent loads, i.e.:

$$ \sum_E (-\mathbf{F}_n^E + \mathbf{f}_n^E + \mathbf{T}_n^E) = 0. \quad (3) $$

### (ii) Identify the weak form of equilibrium of an element as a free body

Element $E$ is subjected to prescribed body forces, and possibly prescribed boundary tractions, and is maintained in a state of equilibrium by the nodal forces $\hat{\mathbf{F}}_n^E$, and $-\mathbf{T}_n^E$, where $\hat{\mathbf{F}}_n^E$ is defined by:

$$ \hat{\mathbf{F}}_n^E = (\mathbf{F}_n^E - \mathbf{f}_n^E). \quad (4) $$

Equilibrium follows from two essential properties of the forces defined in Equation (2), namely that forces $\mathbf{F}_n^E$ for element $E$ are balanced and have a zero resultant, and forces $\mathbf{f}_n^E$, $\mathbf{T}_n^E$ are statically equivalent to the loads applied to element $E$ [14]. The remaining stages of the methodology transform this form of equilibrium to stronger forms. Also, from Equation (3), it is evident that:

$$ \sum_E (-\hat{\mathbf{F}}_n^E - \mathbf{T}_n^E) = 0 \quad (5) $$

so that node $n$ as a free body is in equilibrium under the actions of forces $-\hat{\mathbf{F}}_n^E$, and $\mathbf{T}_n^E$ corresponding to the elements connected to it.

### (iii) Resolve nodal forces at the corners of elements into interactive forces on the sides of elements

The resolving, or "splitting", of corner nodal forces is the key stage to success in this methodology. It is necessary to resolve $\hat{\mathbf{F}}_n^E$, when $n$ refers to a corner node, into components $\hat{\mathbf{P}}_{in}^E$ where $i$ refers to the sides or interfaces of element $E$ which are adjacent to node $n$. For the 2-D membrane there are two sides to consider, and the resolution is not arbitrary, but must satisfy:

$$ \hat{\mathbf{F}}_n^E = \sum_i \hat{\mathbf{P}}_{in}^E \quad \text{and} \quad \sum_E \hat{\mathbf{P}}_{in}^E = 0, \quad (6) $$
where the second summation is over elements E which have node n and interface i in common. The second part of Equation (6) will ensure that interface equilibrium of tractions is also achieved. In the case of a 2-D model this involves only two elements, and solutions to Equation (6) are directly and efficiently provided by a Maxwell force diagram. When n refers to an internal node common to a patch of elements as in Figure 1, then the Maxwell diagram represents Equation (5) without forces $T_n^E$ and is illustrated in Figure 2.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Patch of elements connected to node n.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Maxwell force diagram for node n.</figcaption></figure>

In Figure 2 the vertices correspond to the interfaces in Figure 1 taken in an anticlockwise sequence around node n. This sequence and the insertion of a pole point $P_n$ into the Maxwell diagram enables the resolution to give for example:

$$ \hat{\mathbf{F}}_n^1 = \hat{\mathbf{P}}_{in}^1 + \hat{\mathbf{P}}_{jn}^1; \quad \hat{\mathbf{F}}_n^5 = \hat{\mathbf{P}}_{in}^5 + \hat{\mathbf{P}}_{hn}^5; \quad \text{and} \quad \hat{\mathbf{P}}_{in}^1 + \hat{\mathbf{P}}_{in}^5 = 0. $$

Generally there are two degrees of freedom to this resolution corresponding to the position of the pole point. When node n lies on the boundary Γ various situations occur depending upon the local
boundary conditions [14]. In particular if tractions are specified on the boundary, then the Maxwell diagram includes two vectors $T_n^E$ and the pole point is constrained to be at the vertex formed by the intersection of these two vectors. In this way the resolution of $\hat{\mathbf{F}}_n^E$ complies with the boundary force component $T_n^E$.

(iv) Transform discrete forces on the sides of elements to interactive traction distributions

The replacement of corner forces $\hat{\mathbf{F}}_n^E$ by resolved components $\hat{\mathbf{P}}_{in}^E$ acting at the position of node $n$ but on the adjacent sides does not change the state of equilibrium of an element. In order to define co-diffusive traction distributions $\hat{\mathbf{t}}_i^E$ on side $i$ of element $E$ (i.e. $\sum_{\mathbf{E}} \hat{\mathbf{t}}_i^E = 0$ at all points along a common interface, c.f. the second part of Equation (6)) and maintain element equilibrium, it is now sufficient to replace the discrete forces associated with a side of an element with statically equivalent tractions in a similar way for all elements sharing an interface. There are many ways in which this can be achieved, but the original method [11] is based on satisfying prolongation conditions for conforming elements with a statically admissible stress tensor $\hat{\sigma}^E$ i.e. $\int_{\Omega^E} (\hat{\sigma}^E \cdot \text{grad } \phi_n^E) d\Omega = \int_{\Omega^E} (\sigma_h^E \cdot \text{grad } \phi_n^E) d\Omega$. In effect these conditions establish a unique correspondence between each discrete force and a statically equivalent consistent traction distribution, or mode, using the shape functions $\phi_{in}$ restricted to interface $i$. Then the vector representation $\{\hat{\mathbf{T}}_i^E\}$ of $\hat{\mathbf{t}}_i^E$ using the scalar functions $\phi_{in}$ as a basis for traction distributions is defined by Equation (7).

$$ \{\hat{\mathbf{T}}_i^E\} = [\mathbf{W}]^{-1} \cdot \{\hat{\mathbf{P}}_i^E\}, \qquad (7) $$

where the coefficients of $[\mathbf{W}]$ are defined by $w_{pq} = \int_{\text{side } i} \phi_p^E \phi_q^E d\Gamma$ and indices $p, q$ refer to node numbers on side $i$ of element $E$ assuming a Lagrangian or Serendipity type of element. The components of $\{\hat{\mathbf{P}}_i^E\}$ refer to discrete forces on the nodes of side $i$ resolved into a common direction with the tractions $\hat{\mathbf{t}}_i^E$, e.g. normal or tangential to a straight side. For the 4-noded Lagrange element,

$$ [\mathbf{W}] = \frac{L_i}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}; \quad \text{and} \quad [\mathbf{W}]^{-1} = \frac{2}{L_i} \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix}, $$

and shape functions and traction modes for side $i$ of element $E$ are illustrated in Figure 3.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Linear functions for side displacements and tractions: a) shape functions, b) normal modes of traction</figcaption></figure>
(v) Analyse the subdomain of each element to recover statically admissible stresses

The subdomain of each element now presents a local problem if statically admissible stresses are required as well as equilibrating tractions. Each problem is posed with specified body forces $f^E$ and boundary tractions $\tilde{t}^E$. The internal stress field may, as with the tractions, be defined in many ways depending on the type of replacement element which is chosen to model a subdomain. The original method [11, 12] utilised stress based equilibrium macro-elements of low degree (e.g. 1 or 2). The concept of an equilibrium macro-element free of spurious kinematic modes, and hence capable of transmitting tractions of general degree $p$, has recently been further developed [17] and could be used in the present context. However a more pragmatic approach would allow some relaxation of strict statical admissibility, and this leads to a wider range of replacement elements, e.g. hybrid Trefftz HT-T elements [7] with sufficient degrees of stress fields to obviate the exciting of spurious kinematic modes, and $p$-type displacement elements where degree ($p$ + 3) has been suggested [3] as a pragmatic approach.

It should be noted that this methodology involves a number of non-unique results i.e. the resolution of forces, the definition of codiffusive tractions, and the element stress fields all involve choices, since generally the finite element model is statically indeterminate and even the replacement elements may be hyperstatic. There is therefore scope for optimisation depending on the nature of the problem e.g. error analysis or limit state design.

3. EXTENSION TO 3-D PLATE MODELS WITH NODAL FORCES AND MOMENTS

Plate and/or shell elements are used to model thin "surface" forms of structure. In this section
elements are considered with a single node within the thickness of the model, and surfaces do not
intersect themselves, i.e. a side of an element is common to at most two elements. In this context
of curved surfaces the simplification of using flat elements with straight sides has proved effective
for many problems [19]. Models of plates and shells present two new aspects compared with the
membrane models in 2-D:

(a) The dimension of a nodal force may increase to 3 to represent both in-plane and out-of-plane action. Furthermore the concept of forces and tractions must be extended to include moment vectors and couple distributions, which may also be three dimensional when drilling degrees of freedom are present.

(b) Plate and shell elements are often non-conforming, which makes the use of their shape functions in stages (i) and (iv) problematic.

In principle, however, the five stages of Section 2 still define a general methodology for recovering
equilibrium, provided that nodal force and moment vectors are defined with the essential properties
which imply the equilibrium of nodes and elements as free bodies in stages (i) and (ii). The resolution
stage (iii) now requires separate Maxwell diagrams for both force and moment vectors, and these
diagrams with their pole points may occur in 1-D, 2-D, or 3-D, e.g. transverse shear forces for
flat plates involve 1-D diagrams, moment vectors for plate bending involve 2-D diagrams, and the
general combination of in-plane and out-of-plane actions for folded plates and shells involve 3-D
diagrams. The transformations of discrete forces and moments to distributions of tractions and
couples respectively requires a set of basis functions for an interface. In the absence of conforming
shape functions associated with the elements, use can be made of "imported" shape functions,
e.g. associated with a Lagrange line element mapped onto an interface. These functions should
permit replacement of a discrete vector by statically equivalent distributions. In principle Equation
(7) would be valid for curved interfaces in 3-D, but the integration of the coefficients $w_{pq}$ could
become more complicated to account for a non-linear mapping from a parent straight line element.
Again considerable simplification is introduced here by the use of flat elements having straight
sides. Another complication comes from the replacement of moment vectors by distributions of couples. Distributions of couples about an axis normal to the element surface can be replaced for simplicity by equivalent in-plane tractions in a similar way to the replacement of torsional couples by equivalent Kirchhoff shear forces. In the case when adjacent flat elements exist in different planes, strict codiffusion of traction forces would require couples about both axes normal to an interface to be replaced by equivalent tractions. Thus with reference to Figure 4

$$ \hat{t}_{ix} + \frac{\partial \hat{m}_{ix}}{\partial y} \quad \text{replaces} \quad \hat{t}_{ix}, \quad \hat{t}_{iz} - \frac{\partial \hat{m}_{iz}}{\partial y} \quad \text{replaces} \quad \hat{t}_{iz} \qquad (8) $$

and couple distributions $\hat{m}_{ix}$ and $\hat{m}_{iz}$ are replaced by discrete corner forces equal in value to these distributions at the end points of a side. Finally these discrete forces give rise to additional traction distributions.

<figure><img src="image_6.png" /><figcaption>Fig. 4. Side tractions/couples with reference to local axes for a flat shell element, where $\hat{m}_{ix}$, $\hat{m}_{iy}$, $\hat{m}_{iz}$ refer to torsional, bending and drilling couples per unit length</figcaption></figure>

The recovery of statically admissible stresses in a curved shell subdomain in stage (iv) would also be problematic, but flat subdomains with straight sides are again simpler to analyse. In the case of flat elements it is sufficient to recover statically admissible stress resultants, e.g. moments and shear forces, as for 2-D theories of membranes and plates.

## 4. A PROCEDURE FOR DEALING WITH HIGH DEGREES OF ELEMENT CONNECTIVITY

The resolution of corner nodal forces/moments as expressed by Equation (6) becomes more complex in 3-D models where surface structures intersect, e.g. as in box girders with diaphragms, or where solid elements surround a corner node. In the former case interfaces are connected to or are common with 3 or more shell elements, and in the latter case an element corner node is generally connected to 3 or more 2-D interfaces. Then:

$$ \begin{array}{ll} \text{either} & \displaystyle \sum_{E} \hat{P}_{in}^{E} = 0 \quad \text{for } > 2 \text{ elements } E, \\ \text{or} & \displaystyle \hat{F}_{n}^{E} = \sum_{i} \hat{P}_{in}^{E} \quad \text{for } > 2 \text{ interfaces } i \text{ respectively.} \end{array} \qquad (9) $$
In both cases, although there is no change in principle to the procedure, the use of Maxwell diagrams to obtain solutions to Equation (6) is generally no longer appropriate. Their use depends on resolving nodal forces into only two components which are to act on interfaces connected to two elements. However when the degree of element connectivity is increased solutions to Equation (6) can still be generated in a simple way by exploiting a statical interpretation with graph-theoretic concepts [8, 16]. As with the Maxwell diagram, separating topological and geometrical properties can lead to reduced computational effort when compared to a purely algebraic approach. A graph associated with a node allows direct definition of particular and complementary solutions for the local statical problem. Thus, as with the flexibility method of structural analysis, computational benefits can be realised when the complementary solutions are directly defined based on circuits of a graph.

## 4.1. A reinterpretation of the 2-D case

In Section 2 nodal forces $\hat{\mathbf{F}}_n^E$ have been resolved into two components to satisfy Equation (6) using the concept of a Maxwell diagram. Another interpretation of this resolution allows for a straightforward extension to the case of 3-D models with shell elements or three dimensional elements such as tetrahedra or hexahedra. In 2-D models a patch of elements having an internal node $n$ in common can be considered in a similar way to a frame structure forming a single ring around node $n$ when the elements are structurally connected along their interfaces. If this ring structure is loaded by the set of forces $-\hat{\mathbf{F}}_n^E$, for all $E$, it is in overall equilibrium since these forces maintain node $n$ in equilibrium as a free body (c.f. Equation (5)).

Such a 2-D frame structure is generally three times statically indeterminate. However, in the limit, as the ring is conceptually shrunk around node $n$, the degree of statical indeterminacy reduces to two, and internal forces which form a statical solution are defined by the interface forces $\hat{\mathbf{P}}_{in}^E$ which satisfy Equation (6). The general statical solution for internal forces may be expressed as the combination:

$$ \hat{\mathbf{P}}_{in}^E = {}^o\hat{\mathbf{P}}_{in}^E + {}^c\hat{\mathbf{P}}_{in}^E, \qquad (10) $$

where the first and second terms on the right hand side refer to particular and complementary solutions respectively.

It will be useful, particularly in the 3-D case, to introduce some graph-theoretic concepts which simplify the description of internal force paths around nodes. A directed graph $G$ is defined by two sets of vertices $v^E$ and $v_i$ corresponding to the sets of elements and sides or interfaces respectively, and edges $e_i^E$ corresponding to pairs of vertices $(v^E, v_i)$ when interface $i$ is connected to element $E$. The edges are directed with a positive sense from $v^E$ towards $v_i$. A subgraph $G(n)$ is defined for a corner node $n$ by the vertices, and edge pairs, corresponding to the elements and interfaces connected to node $n$. This is illustrated in Figure 5a where the graph $G(n)$ is embedded in the patch of elements. The vertices of $G$ are denoted by “o” for an “element” vertex $v^E$, and by “x” for an interface vertex $v_i$. For each subgraph $-\hat{\mathbf{F}}_n^E$ is associated with vertex $v^E$, and $\hat{\mathbf{P}}_{in}^E$ is associated with edge $e_i^E$.

### 4.1.1. Particular solution

A particular solution is formed by cutting one of the interfaces, such as $i$, and conducting the internal forces through the remaining interfaces. Force paths thus correspond to paths in the tree of the graph defined by removing one of the edges, say $e_i^1$. Then the internal forces of the particular solution pick up the loads in the order in which they are encountered in the path which starts at vertex $v^1$ and ends at vertex $v_i$. The load force vectors, taken in this order, form the force polygon in Figure 5b. The particular solution is shown in this figure, and this solution is given by Equation (11).
<figure><img src="image_7.png" /><figcaption>Fig. 5. a) Graph G(n) for a patch of elements connected to node n; b) Force polygon with a particular solution; c) Force polygon with a total solution</figcaption></figure>

$$ ^{0}\hat{P}_{jn}^{1} = \hat{F}_{n}^{1} = -^{0}\hat{P}_{jn}^{2}, $$

$$ ^{0}\hat{P}_{kn}^{2} = \hat{F}_{n}^{1} + \hat{F}_{n}^{2} = -^{0}\hat{P}_{kn}^{3}, $$

$$ ^{0}\hat{P}_{mn}^{3} = \hat{F}_{n}^{1} + \hat{F}_{n}^{2} + \hat{F}_{n}^{3} = -^{0}\hat{P}_{mn}^{4}, $$

(11)

$$ ^{0}\hat{P}_{hn}^{4} = \hat{F}_{n}^{1} + \hat{F}_{n}^{2} + \hat{F}_{n}^{3} + \hat{F}_{n}^{4} = -^{0}\hat{P}_{hn}^{5}, $$

$$ ^{0}\hat{P}_{in}^{5} = \hat{F}_{n}^{1} + \hat{F}_{n}^{2} + \hat{F}_{n}^{3} + \hat{F}_{n}^{4} + \hat{F}_{n}^{5} = 0 = -^{0}\hat{P}_{in}^{1}. $$
### 4.1.2. Complementary solutions

Complementary solutions are formed by applying a biaction $\pm \mathbf{R}$ at an interface, say $i$, in the absence of loads $-\hat{\mathbf{F}}_n^E$. The biaction is thus transmitted around a force path corresponding to the single circuit of the graph. Internal forces are defined by Equation (12).

$$ 
\begin{aligned}
& c\hat{\mathbf{P}}_{jn}^1 = c\hat{\mathbf{P}}_{kn}^2 = c\hat{\mathbf{P}}_{mn}^3 = c\hat{\mathbf{P}}_{hn}^4 = c\hat{\mathbf{P}}_{in}^5 = \mathbf{R}; \\
& c\hat{\mathbf{P}}_{in}^1 = c\hat{\mathbf{P}}_{jn}^2 = c\hat{\mathbf{P}}_{kn}^3 = c\hat{\mathbf{P}}_{mn}^4 = c\hat{\mathbf{P}}_{hn}^5 = -\mathbf{R};
\end{aligned}
\qquad (12) $$

### 4.1.3. Total solution

The total solution is formed from the sum in Equation (10) to give for example:

$$ 
\begin{aligned}
\hat{\mathbf{P}}_{jn}^1 &= {}^0\hat{\mathbf{P}}_{jn}^1 + \mathbf{R} = \mathbf{j}\hat{\mathbf{P}}_n; & \hat{\mathbf{P}}_{jn}^2 &= \hat{\mathbf{P}}_n\mathbf{j}, \\
\hat{\mathbf{P}}_{kn}^2 &= {}^0\hat{\mathbf{P}}_{kn}^2 + \mathbf{R} = \mathbf{k}\hat{\mathbf{P}}_n; & \hat{\mathbf{P}}_{kn}^3 &= \hat{\mathbf{P}}_n\mathbf{k}, \\
\hat{\mathbf{P}}_{mn}^3 &= {}^0\hat{\mathbf{P}}_{mn}^3 + \mathbf{R} = \mathbf{m}\hat{\mathbf{P}}_n; & \hat{\mathbf{P}}_{mn}^4 &= \hat{\mathbf{P}}_n\mathbf{m}, \\
\hat{\mathbf{P}}_{hn}^4 &= {}^0\hat{\mathbf{P}}_{hn}^4 + \mathbf{R} = \mathbf{h}\hat{\mathbf{P}}_n; & \hat{\mathbf{P}}_{hn}^5 &= \hat{\mathbf{P}}_n\mathbf{h}, \\
\hat{\mathbf{P}}_{in}^5 &= {}^0\hat{\mathbf{P}}_{in}^5 + \mathbf{R} = \mathbf{i}\hat{\mathbf{P}}_n; & \hat{\mathbf{P}}_{in}^1 &= \hat{\mathbf{P}}_n\mathbf{i}.
\end{aligned}
\qquad (13) $$

where the point $\hat{\mathbf{P}}_n$ in Figure 5c is positioned in the force polygon so that $\mathbf{i}\hat{\mathbf{P}}_n = \mathbf{R} \cdot \hat{\mathbf{P}}_n$ has now the same significance as the pole point in the Maxwell diagram for the 2-D case, and its two degrees of freedom correspond to the two degrees of statical indeterminacy of the ring structure around node $n$.

## 4.2. Extension to the 3-D case

Finite element models in 3-D may contain shell or solid elements, and a set of elements surrounding a corner node forms a three dimensional “ring” structure which itself may contain several rings which allows for more choice in selecting particular and complementary solutions. To simplify presentation the finite element models will be considered as an assembly of shell or solid elements but not a mixture of both. The formation of an associated graph $G$ is now described in a more general way to suit the 3-D case, and it is convenient to use concepts and terms from algebraic topology. The topological properties of a finite element mesh are considered in a cell complex $K$ [15], with 0-, 1-, 2, and 3-cells representing corner nodes, sides or interfaces of shell elements and edges of solid elements, shell elements or interfaces of solid elements, and solid elements respectively. It is assumed that a finite element mesh is regular as regards its connections, so that the cell complex is properly connected. The complex $K$ is now subdivided to $K^1$ with one new 0-cell introduced within each original 1-, 2-, or 3-cell. This subdivision is similar to a mesh refinement, but is only carried out in $K$ for the purposes of defining a suitable graph $G$. Thus for example each 1-cell is subdivided into 2 1-cells, a quadrilateral 2-cell is subdivided into 4 2-cells, and a hexahedral 3-cell is subdivided into 8 3-cells.

In the case of shell models the dimension of $K$ is 2, i.e. no 3-cells are present, and $G$ is defined as the 1-skeleton of $K^1$ less the original 0-cells and their incident 1-cells. This means that the subgraph $G(n)$ associated with node $n$ is the link $Lk(n)$ of 0-cell $n$ in $K^1$ [15]. $Lk(n)$ is usually defined in terms of the neighbourhood $N(n)$ of a 0-cell which consists of all 2-cells incident with the 0-cell together with all their 1- and 0- faces. Then $Lk(n)$ is the 1-skeleton of $N(n)$ less the 0-cell $n$ and the 1-cells incident with it. Two examples of subgraphs $G(n)$ as “links” are illustrated in Figure 6.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Subgraphs $G(n)$ as links $Lk(n)$ of neighbourhoods $N(n)$ of node $n$ for plate/shell models.</figcaption></figure>

In the case of solid models $K$ has dimension 3 and the neighbourhood $N(n)$ in $K^1$ consists of all 3-cells incident with a 0-cell together with all their faces. The link $Lk(n)$ is now the subcomplex consisting of all the 2-cells of $N(n)$ **not** incident with 0-cell $n$ together with all their faces. $|Lk(n)|$ is topologically equivalent to a sphere (for an internal node), or a disc (for a boundary node). For the solid model, the edges of an element are not sites for the application of equilibrating tractions, so that the corresponding 1-cells of $K$, and their subdivisions in $K^1$ are superfluous in defining $G$. Hence $G(n)$ is now defined as the 1-skeleton of $Lk(n)$ less the 0-cells, and their incident 1-cells, which were used in the subdivision of the 1-cells of $K$. Two examples of $Lk(n)$ and $G(n)$ for patches of hexahedral and tetrahedral elements are illustrated in Figures 7 and 10. $G$ becomes the union of $G(n)$ for all $n$.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Patch of hexahedral elements showing partial views of $Lk(n)$ and $G(n)$</figcaption></figure>

It should be noted that nodes on the edges of **solid** elements do present sites where the resolution of nodal forces is required since these forces must be transferred to 2-D interfaces. However these
nodes can be treated in stage (iii) in exactly the same way as described for membrane corner nodes in Section 4.1 except that the forces are 3-dimensional not 2-dimensional. In this case the subgraph $G(n)$ is now defined as the link $Lk(n)$ when the neighbourhood $N(n)$ is restricted to the subcomplex consisting of all 2-cells which are incident with 0-cell $n$ and the 0-cells introduced within the original 3-cells of $K$. Hence for an internal edge node, $G(n)$ forms one circuit as in Figure 5a, and the concept of a 3-D Maxwell diagram becomes appropriate again.

The subgraph $G(n)$ defines the components of the local "ring structure" for node $n$, with vertex $v^E$ representing element $E$ and load $-\hat{\mathbf{F}}_n^E$, vertex $v_i$ representing interface $i$ where structural interaction between elements occurs, and edge $e_i^E$ represents the interactive force $\hat{\mathbf{P}}_{in}^E$. A statical solution for the ring structure is again expressed in terms of particular and complementary solutions, both of which involve the cyclomatic number $\mu(n)$ of $G(n)$ defined by Equation (14).

$$ \mu(n) = (n_e - n_v + 1), \qquad (14) $$

where $n_e$ and $n_v$ are the numbers of edges and vertices respectively in $G(n)$.

### 4.2.1. Particular solution

A particular solution is formed by defining force paths to transmit forces such as $-\hat{\mathbf{F}}_n^E$ to a reference element $B$ say. A path from $v^E$ to $v^B$ is defined in a tree of $G(n)$, a tree being selected by removing $\mu(n)$ edges so that the vertices remain connected. The incidence of an edge in such a path is recorded in the $(n_v^E - 1) \times n_e$ incidence matrix $[\mathbf{T}]$, where $n_v^E$ is the number of vertices representing elements. Then a particular statical solution can be written in the form of Equation (15).

$$ \{^{o}\hat{\mathbf{P}}_n\} = [\mathbf{B}_o]\{-\hat{\mathbf{F}}_n\}, \qquad (15) $$

where $\{^{o}\hat{\mathbf{P}}_n\}$ contains $3n_e$ components of interface forces, $\{-\hat{\mathbf{F}}_n\}$ contains $3(n_v^E - 1)$ components of element forces which excludes the force on element $B$, and $[\mathbf{B}_o]$ is pattern equivalent to $[\mathbf{T}]^T$, i.e. each coefficient $t_{ij} = \pm 1$ or 0 is replaced by $t_{ij} [\mathbf{I}]$ where $[\mathbf{I}]$ is the $3 \times 3$ unit matrix. $[\mathbf{B}_o]^T$ may thus be expressed as the Kronecker product $[\mathbf{T}] \otimes [\mathbf{I}]$.

### 4.2.2. Complementary solutions

Complementary solutions are formed by transmitting biactions around force paths corresponding to $\mu(n)$ independent circuits from the graph $G(n)$. The incidence of an edge in a circuit is recorded in the $\mu(n) \times n_e$ incidence matrix $[\mathbf{C}]$, and the complementary solution appears in the form of Equation (16).

$$ \{^{c}\hat{\mathbf{P}}_n\} = [\mathbf{B}]\{\mathbf{R}_n\}, \qquad (16) $$

where $\{^{c}\hat{\mathbf{P}}_n\}$ contains $3n_e$ components of interface forces, $\{\mathbf{R}_n\}$ contains $3\mu(n)$ components of biactions, and $[\mathbf{B}]$ is pattern equivalent to $[\mathbf{C}]^T$ in a similar way to the pattern equivalence of $[\mathbf{B}_o]$ and $[\mathbf{T}]^T$. $[\mathbf{B}]^T$ may thus also be expressed as the Kronecker product $[\mathbf{C}] \otimes [\mathbf{I}]$.

### 4.2.3. Total solution

The total solution has the form:

$$ \{\hat{\mathbf{P}}_n\} = \{^{o}\hat{\mathbf{P}}_n\} + \{^{c}\hat{\mathbf{P}}_n\} = [\mathbf{B}_o]\{-\hat{\mathbf{F}}_n\} + [\mathbf{B}]\{\mathbf{R}_n\}. \qquad (17) $$

This solution has $3\mu(n)$ variables as components of $\{\mathbf{R}_n\}$ whose values can be freely chosen. However it may be desirable for the total solution to be optimised in the sense of being close to
target forces $\{\bar{\mathbf{P}}_n\}$. For the solid model the target forces may be taken as the mean of the discrete interface forces for two adjacent elements, based on forces consistent with the tractions directly equilibrating with the stresses in the displacement elements. This approach has been advocated for 2-D membrane elements [12]. The statical solution is brought close to the target forces by minimising $\{\bar{\mathbf{P}}_n - \hat{\mathbf{P}}_n\}^T \{\bar{\mathbf{P}}_n - \hat{\mathbf{P}}_n\}$ with respect to the variables, i.e. by solving the Equation (18):

$$[\mathbf{B}]^T [\mathbf{B}] \{\mathbf{R}_n\} = [\mathbf{B}]^T \{\bar{\mathbf{P}}_n - {}^o\hat{\mathbf{P}}_n\}. \qquad (18)$$

In these equations there is no coupling between the three components of the forces when referred to say orthogonal Cartesian axes, so that three independent sets of equations exist in the form:

$$[\mathbf{C}] [\mathbf{C}]^T \{R_{nx}\} = [\mathbf{C}] \{\bar{P}_{nx} - {}^o\hat{P}_{nx}\} \qquad (19) \\ \mu(n) \times \mu(n)$$

where “x” refers to a typical component axis. The matrix $[\mathbf{C}] [\mathbf{C}]^T$ is symmetric, positive definite, low in dimension, and is generally sparse. Furthermore, since it only depends on the local topological properties of a patch of elements, it may be repeated for many nodes of a uniform mesh. Thus, as pointed out by Ainsworth [1] in another method related to the recovery of self equilibrated fluxes and involving local topological matrices, considerable savings in computational effort should be possible from solving such sets of equations.

A very low computational effort should be necessary for selecting trees and independent circuits from the local subgraphs $G(n)$. A number of efficient algorithms exist for tree selection [4], and independent circuits can always be selected as the fundamental circuits associated with a tree. However each subgraph is planar and other algorithms for the selection of independent circuits become appealing. In the case of a solid model, $G(n)$ is embedded into a sphere or disc in the form of $|Lk(n)|$ in the process of its formation. This embedding can be made use of since it divides the sphere or disc into a number of “regional” discs which are bounded by “regional” circuits [5]. When the embedding is into a sphere, the number of regional circuits is $(\mu(n) + 1)$ and the removal of any one circuit leaves the correct number of independent circuits. The two examples of solid models illustrated in Figures 7 and 10 lead to the subgraphs $G(n)$ shown in Figures 8 and 11 respectively. When the embedding is into a disc, the regional circuits number $\mu(n)$ and they are already independent. In the case of a shell model, $G(n)$ as for example in Figure 6 can be embedded into a sphere $S(n)$ which is the boundary of a solid ball formed by expanding the 0-cell, or node $n$.

<figure><img src="image_10.png" /><figcaption>Fig. 8. Subgraph $G(n)$ for a patch of 8 hexahedral elements, with independent regional circuits and a circuit emphasized as a Hamiltonian circuit when it includes the “element” vertices only</figcaption></figure>
### 4.2.4. Examples

The two examples of solid models are now discussed further, and they serve to illustrate the potential benefits from the use of a graph-theoretic approach. In such models the numbers of variables becomes significant for even "small" problems, so that the quest for efficient algorithms has a high priority. The two examples concern patches of hexahedral and tetrahedral elements connected to an internal corner node $n$.

The corner node of each element has three adjacent interfaces, and each interface is connected to two adjacent elements. Denoting the number of vertices in $G(n)$ corresponding to interfaces by $n_v^i$, this means that the characteristics of $G(n)$ are as follows:

$$n_e = 3n_v^E = 2n_v^i \Rightarrow \mu(n) = n_e - n_v + 1 = \left(\frac{1}{2}n_v^E + 1\right). \qquad (20)$$

In the context of the minimisation problem solved by Equations (18) and (19), it is of interest to compare $\mu(n)$ with $\lambda(n)$ where $\lambda(n)$ is the number of Lagrange multipliers in an earlier alternative formulation of the problem [13]. In that formulation there are $n_v^i$ statical variables $\hat{P}_{in}$ representing interface biactions, i.e. $\hat{P}_{in}^E = \hat{P}_{in} = -\hat{P}_{in}^M$ when elements $E$ and $M$ share interface $i$. The minimisation problem is equivalent to minimising $\sum_i (\bar{P}_{in} - \hat{P}_{in})^2$ subject to $\lambda(n)$ independent equilibrium constraints on the force variables $\hat{P}_{in}$. This leads to a total of $(n_v^i + \lambda(n))$ independent vector equations. Eliminating the statical variables leaves :

$$\lambda(n) = (n_v^E - 1) \qquad (21)$$

vector equations for the Lagrange multipliers. Thus the approach based on circuits reduces the number of simultaneous equations by almost one half.

#### (a) Hexahedral elements

Figure 7 illustrates a patch of 8 hexahedral elements, so that $G(n)$ has the following characteristics:

$$n_v^E = 8; \quad n_v^i = 12; \quad n_v = 20; \quad n_e = 24; \quad \mu(n) = 24 - 20 + 1 = 5.$$

Thus for each component of $\hat{P}_{in}^E$ there are 5 independent statical variables, and these can be associated with the 5 independent circuits of $G(n)$ shown in Figure 8. The minimisation problem solved by Equation (19) involves the 5 × 5 symmetric matrix:

$$[C][C]^T = 2 \begin{bmatrix} 4 & -1 & 0 & -1 & -1 \\ -1 & 4 & -1 & 0 & -1 \\ 0 & -1 & 4 & -1 & -1 \\ -1 & 0 & -1 & 4 & -1 \\ -1 & -1 & -1 & -1 & 4 \end{bmatrix}$$

which is the same for every internal corner node, and thus its inverse or factorisation is only required once. The alternative approach involves 7 Lagrange multipliers. Thus the use of circuits gives a small reduction in the number of simultaneous equations to be solved. The second example shows a much bigger reduction.

#### (b) Tetrahedral elements

It may be more appropriate to split or decompose a hexahedral domain into 5 or 6 tetrahedral elements [19]. A decomposition into 5 is shown in Figure 9, and the numbers of tetrahedral elements connected to the 8 corner nodes is indicated as 1 or 4. Thus if the model in example (a) is
<figure><img src="image_11.png" /><figcaption>Fig. 9. Decomposition of a hexahedron into 5 tetrahedrons</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Patch of tetrahedral elements showing partial views of $Lk(n)$ and $G(n)$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Subgraph $G(n)$ for a patch of 32 tetrahedral elements connected to node $n$</figcaption></figure>
remeshed with tetrahedral elements, then the local topology at the internal node $n$ depends on the orientations of the decompositions. If the orientations are all the same, then:

$$n_v^E = 4 \times 1 + 4 \times 4 = 20 \Rightarrow \mu(n) = 11, \lambda(n) = 19.$$

The maximum number of connections to node $n$ occurs when the orientations are different, and this is illustrated in Figure 10 with the corresponding subgraph $G(n)$ shown in Figure 11. Then:

$$n_v^E = 8 \times 4 = 32 \Rightarrow \mu(n) = 17, \lambda(n) = 31.$$

The potential reduction shown in the number of variables should thus lead to a significant improvement in computational efficiency of the traction recovery procedure particularly for solid models.

## Hamiltonian circuits

Although 3-D models with higher degrees of element connectivity no longer involve the concept of a Maxwell diagram in general, an interesting property exists when a graph is Hamiltonian: that is when a simple circuit passes through all its vertices. In this case a tree can be selected by removing one edge from a Hamiltonian circuit, and all the vertices of the tree are connected to two edges at most. This is similar to the 2-D case, and means that a particular solution could be based on a single 3-D Maxwell diagram, and it would then have the simplicity of construction indicated in Figure 5b and Equation (11). The previous example (a) of the hexahedral elements has a Hamiltonian graph, this property is demonstrated in Figure 8. Another example of a Hamiltonian graph $G(n)$ is shown in Figure 12. Unfortunately it appears that at present no simple general test exists to determine whether even a planar graph is Hamiltonian or not [4], so this remains an open question for the second example with tetrahedral elements illustrated in Figures 10 and 11!

<figure><img src="image_14.png" /><figcaption>Fig. 12. Subgraph $G(n)$ embedded: a) on the edges of a dodecahedron, and b) in the plane with a circuit emphasized as a Hamiltonian circuit when it includes the “element” vertices only</figcaption></figure>

## 5. CONCLUSIONS

* the methodology and procedures described and illustrated in this paper enable strong equilibrium to be recovered from 2-D and 3-D displacement models exhibiting only weak equilibrium,

* only small localised problems are involved, no global reanalysis is necessary,
* the use of graph-theoretic concepts has the potential to simplify the implementation of traction recovery and thereby to increase the computational efficiency of the methodology,

* further work is required in this presentation to cover irregular meshes for plate and solid models,

* further work is required to address the associated problems of selecting the optimum equilibrating solutions, particularly for plate and solid models. The question of optimisation depends on the problem being considered, and the criteria may be different for error analysis, limit state design etc.

## ACKNOWLEDGEMENTS

The authors wish to express their thanks for support from C.N.R.S. and the British Council Alliance project PN95.157 which made collaboration possible.

## REFERENCES

[1] M. Ainsworth, J.T. Oden. A posteriori error estimators for second order elliptic systems, Part 2: an optimal order process for calculating self equilibrating fluxes. *Computers and Mathematics*, **26**: 75-88, 1993.

[2] M. Ainsworth, J.T Oden, W. Wu. A posteriori error estimation for *h*-p approximations in elastostatics. *Applied Numerical Mathematics*, **14**: 23-54, 1994.

[3] I. Babuska, T. Strouboulis, A. Mathur, C.S. Upadhyay. Pollution error in the h-version of F.E.M. and the local quality of a posteriori error estimators. Technical Note BN1163. University of Maryland, 1994.

[4] A. Gibbons. *Algorithmic graph theory*. Cambridge University Press, 1985.

[5] J.C.de C. Henderson, E.A.W. Maunder. A problem in applied topology: on the selection of cycles for the flexibility analysis of skeletal structures. *J. Inst. Maths. Applic.*, **5**: 254-269, 1969.

[6] B. Jacobsen. The Slepner accident and its causes. In: J. Robinson, ed., *FEM TODAY and the FUTURE*, 102-108. Robinson and Associates, 1993.

[7] J. Jirousek, A.P. Zielinski. Dual hybrid-Trefftz element formulation based on independent boundary traction frame. *Int. J. Num. Meth. in Eng.*, **36**: 2955-2980, 1993.

[8] A. Kaveh, *Structural Mechanics: Graph and Matrix Methods*, Research Studies Press Ltd., 1992.

[9] D.W. Kelly, J.D. Isles. Procedures for residual equilibration and local error estimation in the finite element method. *Comm. Appl. Numer. Meth.*, **5**: 497-505, 1989.

[10] T. Kvamsdal, K.M. Mathisen, K.M. Okstad, H.B. Aas. Techniques for reliable calculation of sectional forces. In: J. Robinson, ed., *FEM TODAY and the FUTURE*, 450-456. Robinson and Associates, 1993.

[11] P. Ladeveze. Nouvelle procédure d'estimation d'erreur relative à la méthode des éléments finis et d'application. *Actes des Journées Eléments Finis*, Rennes, 1977.

[12] P. Ladeveze, J.P Pelle, Ph. Rougeot. Error estimation and mesh optimization for classical finite elements. *Eng. Comp.*, **8**: 69-80, 1991.

[13] P. Ladeveze, P. Marin, J.P Pelle, J.L. Gastine. Accuracy and optimal meshes in finite element computation for nearly incompressible materials. *Comp. Meth. in Appl. Mech. and Eng.*, **94**: 303-315, 1992.

[14] P. Ladeveze, E.A.W. Maunder. A general method for recovering equilibrating element tractions. *Comp. Meth. in Appl. Mech. and Eng.*, **137**: 111-151, 1996.

[15] C.R.F. Maunder, *Algebraic Topology*, Van Nostrand Reinhold, 1970.

[16] E.A.W. Maunder, G.J. Savage. A graph-theoretic model for finite elements with side variables. *Civil Eng. Syst.*, **11**: 111-141, 1994.

[17] E.A.W. Maunder, J.P. Moitinho de Almeida, A.C.A. Ramsay. A general formulation of equilibrium macro-elements with control of spurious kinematic modes. *Int. J. Num. Meth. in Eng.*, **39**: 3175-3194, 1996.

[18] H. Ohtsubo, M. Kitamura. Element by element a posteriori error estimation and improvement of stress solutions for two-dimensional elastic problems. *Int. J. Num. Meth. in Eng.*, **29**: 223-244, 1990.

[19] O.C. Zienkiewicz, R.L. Taylor, *The Finite Element Method*, 4th Ed., Vols 1 and 2. McGraw-Hill, 1989 and 1991.