Grid Discretization

Among three-dimensional constant strain-rate elements, tetrahedra have the advantage of not generating hourglass deformations (i.e., deformation patterns created by combinations of nodal velocities producing no strain rate, and thus no nodal force increments). However, when used in the framework of plasticity, these elements do not provide for enough modes of deformation (see Nagtegaal et al. 1974). In particular situations, for example, they cannot deform individually without change of volume as required by certain important constitutive laws. In those cases, the elements are known to exhibit an overly stiff response compared to what is expected from theory. To overcome this problem, a process of mixed discretization is applied in FLAC3D.

The principle of the mixed discretization technique is to give the element more volumetric flexibility by proper adjustment of the first invariant of the tetrahedra strain-rate tensor. [*] Application of the technique differs depending on whether primary discretization of the body is done into hexahedral zones or into tetrahedra. In the approach pioneered by Marti and Cundall (1982), a coarser discretization in zones is superposed to the tetrahedral discretization, and the first strain-rate invariant of a particular tetrahedron in a zone is evaluated as the volumetric-average value over all tetrahedra in the zone. The method is illustrated in Figure 1. In the particular mode of deformation sketched there, individual constant strain-rate elements will experience a volume change incompatible with a theory of incompressible plastic flow. In this example, however, the volume of the assembly of tetrahedra (i.e., the zone) remains constant, and application of the mixed discretization process allows each individual tetrahedron to reflect this property of the zone, hence reconciling its behavior with that predicted by the theory.

Mixed Discretization for a Hexahedral Grid

Figure 1: Mode of deformation for which mixed discretization would be most efficient.

In FLAC3D, a zone corresponds to an assembly of \(n_t\) tetrahedra, as illustrated in Figure 2 for the case \(n_t\) = 5. Consider a particular zone where the strain-rate tensor of a tetrahedron locally labeled \(l\) in that zone is first estimated, and then decomposed into deviatoric and volumetric parts:

(1)\[ {\xi}_{ij}^{[l]} = {\eta}_{ij}^{[l]} + {{{\xi}^{[l]}}\over{3}}{\delta}_{ij}\]

where \(η^{[l]}\) is the deviatoric strain-rate tensor, and \(ξ^{[l]}\) is the strain-rate first invariant,

(2)\[ {\xi}^{[l]} = {\xi}_{ii}^{[l]}\]

Figure 2: An 8-node zone with 2 overlays of 5 tetrahedra in each overlay.

The first invariant for the zone is then calculated as the volumetric average value of the first invariant over all tetrahedra in the zone:

(3)\[ {\xi}^z = {{\sum_{k=1}^{n_t} {\xi}^{[k]} V^{[k]}}\over{\sum_{k=1}^{n_t}V^{[k]}}}\]

where \(V^{[k]}\) is the volume of tetrahedron \(k\). F inally, the tetrahedron strain-rate tensor components are calculated from

(4)\[ {\xi}_{ij}^{[l]} = {\eta}_{ij}^{[l]} + {{{\xi}^z}\over{3}} {\delta}_{ij}\]

Dilatant constitutive laws will produce changes in mean normal stress when yielding occurs. For a consistent technique, the first invariant of the stress tensor, derived after application of the strain-rate increment, must also be evaluated as a volumetric average for the zone. In this process, the stress tensor of a particular tetrahedron \(l\) in a zone is first estimated and decomposed into deviatoric and volumetric parts:

(5)\[ {\sigma}_{ij}^{[l]} = s_{ij}^{[l]} + {\sigma}^{[l]} {\delta}_{ij}\]

where [\(s\)]\(^{[l]}\) is the deviatoric strain-rate tensor, and \(σ^{[l]}\) is the mean normal stress.

(6)\[ {\sigma}^{[l]} = {{1}\over{3}} {\sigma}_{ii}^{[l]}\]

The first invariant for the zone is calculated as the volumetric average value over all tetrahedra in the zone:

(7)\[ {\sigma}^z = {{\sum_{k=1}^{n_t}{\sigma}^{[k]} V^{[k]}}\over{\sum_{k=1}^{n_t}V^{[k]}}}\]

Finally, the tetrahedron stress-rate tensor components are calculated using

(8)\[ {\sigma}_{ij}^{[l]} = s_{ij}^{[l]} + {\sigma}^z {\delta}_{ij}\]

In FLAC3D, the discretization process starts with the coarser grid: zones are defined and then discretized (internally) into tetrahedra. An eight-noded zone, for instance, can be discretized into two (and only two) different configurations of five tetrahedra (corresponding to overlay 1 and 2 in Figure 2. The calculation of nodal forces (based on evaluation of strain rates and stresses) can be carried out using one overlay or a combination of two overlays. The advantage of the two-overlay approach is to ensure symmetric zone response for symmetric loading. In this case, mixed discretization is carried out over the combination of two overlays, and nodal forces computations are evaluated by averaging over the two overlays.

Nodal Mixed Discretization for a Tetrahedral Grid

Nodal mixed discretization (NMD) is a variation on the mixed discretization scheme, in which averaging of the volumetric behavior is carried out on a node basis instead of a zone basis. The procedure is applied on the triangle or tetrahedral-based mesh; it does not require the assembly of elements into zones. Also, the constitutive model is called on an element basis, as usual (no call is made on a node basis for the volumetric behavior, as in the average nodal pressure (ANP) formulation described by Bonet and Burton 1998). The procedure involves a nodal mixed discretization on stress. Both steps are considered below.

First we recall the general calculation sequence embodied in FLAC and FLAC3D:

1. Nodal forces are calculated from stresses, applied loads and body forces (velocity and displacement vary linearly; stress and strain are constant within an element).
2. The equations of motion are invoked to derive new nodal velocities and displacements.
3. Element strain rates are derived from nodal velocities.
4. New stresses are derived from strain rates, using the material constitutive law.

In the NMD technique, the calculation sequence is respected. However, an averaging procedure is carried out on stresses (end of step 4), as described below.

First, nodal values for the volumetric stress \(\sigma\) are calculated as the weighted average of the surrounding element values:

(9)\[ \sigma_n = {{\sum_{e = 1}^{m_n} \sigma V_e} \over {\sum_{e = 1}^{m_n} V_e}}\]

where \(m_n\) are the elements surrounding node \(n\), and \(V_e\) is the volume of element \(e\). After the nodal values \(\sigma_n\) are obtained, a mean value for the element \(\bar{\sigma}\) is calculated by taking the average of nodal values:

(10)\[ \bar{\sigma} = {1 \over d} \sum_{n = 1}^d \sigma_n\]

where \(d\) = 3 for a triangle in FLAC and 4 for a tetrahedra in FLAC3D.

Finally, the stresses calculated by the constitutive model are corrected by substituting \(\bar{\sigma}\) for \(\sigma\) in the returned value:

(11)\[ \sigma_{ij} => \sigma_{ij} + (\sigma - \bar{\sigma}) \delta_{ij}\]

It is important to note that the averaging manipulations involved in the NMD technique are done independent of the constitutive model formulation, which remains unaffected.

Also, because of the way the averaging procedure is carried out, the resulting tetrahedral formulation gives the correct behavior in patch test situations where \(\sigma\) is uniform, and thus equal, for all elements.

The nodal averaging procedure involved in the NMD technique removes the excessively constrained kinematics of the linear velocity element.

[*]	This invariant gives a measure of the rate of dilation of the constant strain-rate tetrahedron.