Numerical Implementation

Body Discretization

By default, discretization of the volume under study is done into hexahedral zones. The user can also import a tetrahedral mesh (generated by an external mesher) into FLAC3D. The two types of body discretization are addressed:

Hexahedral Meshing

In FLAC3D, the general discretization of the body into hexahedral zones is performed by the user. Several available procedures to help in this task are described in defcrossrefFilename{ftd103}refvolI. Each zone is discretized automatically by the code into sets of tetrahedra. The user can decide to carry out the calculation using one overlay or a combination of two overlays. The use of two overlays is recommended in regions of the medium where high gradients of stresses and deformations are expected. By default, two overlays are prescribed for all zones in a model. Also, the technique of mixed discretization is applied to overcome the overly stiff behavior of uniform strain tetrahedra during plastic flow.

Tetrahedral Meshing

The NMD calculation carried out on a tetrahedral grid imported into FLAC3D uses only one overlay. Implementation of the nodal mixed discretization technique on strain is done using three calculation loops, performed just before calling the constitutive model:

1. One loop on zones to distribute volumetric strain rate (as calculated from gridpoint velocities) and area/volume to gridpoints.
2. A loop on gridpoints to evaluate the volume-weighted average of volumetric strain rate at the gridpoints.
3. A loop on zones to calculate an average of volumetric strain rate using the gridpoint values, and assign the value to the zone.

Note that advantage can be taken of existing calculation loops to perform the task described above.

In case constitutive models which rely on an elastoplastic volumetric behavior are used, similar loops (based on stress) are used, after calling the constitutive model, to implement the nodal mixed discretization on stress.

Initial and Boundary Conditions

The boundary conditions of the problem consist of surface tractions, concentrated loads and displacements. In addition, body forces may be given and initial stress conditions imposed. For implementation in the code, all stresses and nodal velocities are initially set to zero; then, initial stresses are specified. Concentrated loads are specified at given surface nodes, and imposed boundary displacements are prescribed in terms of nodal velocities. Body forces and surface tractions are transformed internally into a set of statically equivalent nodal forces. This constitutes the initial state of the numerical scheme.

Main Calculation Steps

FLAC3D uses an explicit “time-marching” finite difference solution scheme; for every timestep, the calculation sequence can be summarized:

1. New strain rates are derived from nodal velocities.
2. Constitutive equations are used to calculate new stresses from the strain rates and stresses at the previous time.
3. The equations of motion are invoked to derive new nodal velocities and displacements from stresses and forces.

The sequence is repeated at every timestep, and the maximum out-of-balance force in the model is monitored. This force will either approach zero, indicating that the system is reaching an equilibrium state, or it will approach a constant, nonzero value, indicating that a portion (or all) of the system is at steady-state (plastic) flow of material. The calculation may be interrupted at any point, in order to analyze the solution.

A summary of equations (derived in ref) that are used in the calculation sequence is presented below, for clarity. In the notation, superscripts are used to particularize a node or tetrahedron only when needed in the local context. Also, as before, the notation [[f]]<l> represents the sum of contributions of type f at global node l of all tetrahedra having that node in common.

Strain-Rate Calculation

Starting from a known velocity field, the components of the strain-rate tensor are computed, for each tetrahedron in a zone, using the finite difference formulation Eqref{109E1015}:

(1)\[ {\xi}_{ij} = - {{1}\over{6V}} \sum_{l = 1}^{4} \left( v_i^l\ n_j^{(l)} + v_j^l\ n_i^{(l)}\right) S^{(l)}\]

For the default hexahedral grid, a procedure of mixed discretization is then applied, and new diagonal strain-rate tensor components are calculated from Eqrefs{109E1071} and eqref{109E1072}. For a particular tetrahedron, now locally labeled l in the zone, these equations yield

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

where ξz, the average value for the zone of the first strain-rate invariant, is calculated from

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

where nt is the total number of tetrahedra involved in the zone computation, and V[k] is the volume of a tetrahedron labeled k.

For a tetrahedral grid, a technique of nodal mixed discretization on strain is applied.

Stress Calculation

The constitutive equations are used in their incremental form H*ij to calculate stress increments for each tetrahedron in a zone (see Eqrefs{109E1045}, eqref{109E1046} and eqref{109E1047}):

(4)\[ \Delta {\sigma}_{ij} = \Delta {\check{\sigma}}_{ij} + \Delta {\sigma}_{ij}^C\]

where

(5)\[ \Delta {[\check \sigma}_{ij}] = H_{ij}^* \ ({\sigma}_{ij}, \Delta _{ij})\]

and \(\epsilon_{ij}\) is computed (using Eqref{109E1015}) from

(6)\[ \Delta _{ij} = - {{\Delta t}\over{6V}} \sum_{l = 1}^{4}\left( v_i^l\ n_j^{(l)} + v_j^l\ n_i^{(l)}\right) S^{(l)}\]

The stress correction ΔσCij is not taken into consideration in small-strain mode. In large-strain mode, its value is given by Eqrefs{109E1048} and eqref{109E1049}:

(7)\[ \Delta {\sigma}_{ij}^C = ({\omega}_{ik}\ {\sigma}_{kj} - {\sigma}_{ik}\ {\omega}_{kj}) \Delta t\]

where

(8)\[ {\omega}_{ij} = - {{1}\over{6V}} \sum_{l=1}^4 \left( v_i^l\ n_j^{(l)}- v_j^l\ n_i^{(l)} \right) S^{(l)}\]

New stress values are then derived by addition of the stress increments.

For the default hexahedral grid, a technique of mixed discretization is applied, and the diagonal components of the stress-tensor components are adjusted using Eqrefs{109E1075} and eqref{109E1076}. For a particular tetrahedron, locally labeled l in the zone, these equations are

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

and

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

For a tetrahedral grid, a procedure of nodal mixed discretization on stress is used.

Nodal Mass Calculation

The tetrahedron contribution ml to the nodal mass at local node l is evaluated using Eqref{109E1063}:

(11)\[ m^l ={{{\alpha}_1}\over{9V}} max \left( {\left[ n_i^{(l)} S^{(l)} \right] }^2 {, i=1,3}\right)\]

noting that the magnitude of the timestep is set equal to unity.

The nodal mass M<l> at global node l is computed as the sum of contributions Equation (11) at the node of all tetrahedra having that node in common:

(12)\[ M^{<l>} = {[[m]]}^{<l>}\]

In the small-strain mode, the nodal mass calculation is only carried out once before cycling. In large-strain mode, values are updated every 10 steps.

Out-of-Balance Force Calculation

The tetrahedron contribution [p]l to the out-of-balance force at local node l is computed from the stresses and body forces using

(13)\[ p_i^l = {{1}\over{3}} {\sigma}_{ij}\ n_j^{(l)}\ S^{(l)} + {{1}\over{4}} \rho\ b_i\ V\]

(see Eqrefs{109E1021} and eqref{109E1040}).

The sum of contributions Equation (13) at global node l of all tetrahedra having that node in common is evaluated. Averaging over two overlays is carried out when appropriate, and the contributions [P]<l> of applied loads and concentrated forces are added, thus providing the expression for the out-of-balance force [F]<l> at global node l:

(14)\[ F_i^{<l>} = {[[ p_i ]]}^{<l>} + P_i^{<l>}\]

The out-of-balance force is monitored to detect whether the system has reached a state of equilibrium or steady flow. The damping term is evaluated from Eqref{109E1065} as

(15)\[\mathcal{F}_{(i)}^{<l>} = -\alpha \left| F_i^{<l>} \right| {\rm sign} \left( v_{(i)}^{<l>} \right)\]

and is added to the out-of-balance force.

Velocity and Displacement Calculations

Eqref{109E1042} is invoked, taking the damping into account, to derive new nodal velocities from known out-of-balance forces:

(16)\[v_i^{<l>}(t + {\Delta t \over 2}) = v_i^{<l>}(t - {\Delta t \over 2}) + {{\Delta t}\over{{M}^{<l>}}} \left( F_i^{<l>} + \mathcal{F}_i^{<l>} \right)\]

In turn, nodal displacements are evaluated using Eqref{109E1044}:

(17)\[ u_i^{<l>}(t + \Delta t) = u_i^{<l>}(t) + \Delta t v_i^{<l>} (t + {\Delta t \over 2})\]

Geometry Update Calculation

In large-strain mode, node locations are updated using Eqref{109E1043}, and the grid then deforms accordingly. For global node l, we have

(18)\[ x_i^{<l>}(t + \Delta t) = x_i^{<l>}(t) + \Delta t v_i^{<l>} (t + {\Delta t \over 2})\]

In small-strain mode, the geometry adjustment is neglected.

Energy Calculation in FLAC3D

Elastic strain energy and dissipated plastic energy can be tracked for zones containing a mechanical model. FLAC3D uses an incremental solution procedure: the equations of motion (at the gridpoints) and the stress-strain calculations (at the zones) are solved every timestep. In the stress-strain equations, the incremental change in energy components is determined and accumulated as the system attempts to reach equilibrium.

Energy Dissipation in Zones through Plastic Work

Several plasticity models are available in FLAC3D which can describe the deformability of the zones. Energy is dissipated through plastic work as the zones undergo irreversible deformation. The strain in any zone can be divided into elastic and plastic parts.

The elastic strain energy (decomposed into deviatoric and volumetric components) is given by

(19)\[ W_e = {V \over 2} \left({{\sigma_{ij}^d \sigma_{ij}^d} \over 2G} +{{\bar \sigma \bar \sigma} \over K}\right)\]

where V is sub-tet volume, σdij is deviatoric stress and \(\bar \sigma\) is mean stress.

The change in elastic strain energy before and after the stress-strain equations are solved in a constitutive model is given by

(20)\[ W_e = W_e - W_e^{\prime}\]

Where We is the current elastic strain energy and Wʹe is the previous elastic strain energy. Equation (20) can be decomposed into shear and volumetric components: ΔWes and ΔWev.

The total shear energy change can be found from the total deviatoric strain and deviatoric stress:

(21)\[\Delta W_{Ts} = {V \over 2} \left[ \left(\sigma_{11} + \sigma^{\prime}_{11}\right)e_{11} + \left(\sigma_{22} + \sigma^{\prime}_{22}\right)e_{22} + \left(\sigma_{33} + \sigma^{\prime}_{33}\right)e_{33} + 2\left(\sigma_{12} + \sigma^{\prime}_{12}\right)e_{12} + 2\left(\sigma_{13} + \sigma^{\prime}_{13}\right)e_{13} + 2\left(\sigma_{23} + \sigma^{\prime}_{23}\right)e_{23} \right]\]

The total volumetric change can be found from the total mean strain and mean stress:

(22)\[ \Delta W_{Tv} = {3V \over 2} \left(\bar \sigma + \bar \sigma^{\prime}\right)\bar e\]

The total shear plastic work dissipated during a timestep is the difference between the total shear energy and the elastic shear energy change at any timestep:

(23)\[ \Delta W_{ps} = \Delta W_{Ts} - \Delta W_{es}\]

Similarly, the total volumetric plastic work dissipated during a timestep is

(24)\[ \Delta W_{pv} = \Delta W_{Tv} - \Delta W_{ev}\]

Each zone contains overlays of constant stress, strain tetrahedra. A volumetric average of the above energy components is calculated for the zone for an average energy dissipation increment for that zone. Each zone has four additional variables to accumulate dissipated plastic energy: Wpv, Wps, Wev, Wes.

Energy calculations in FLAC3D can be accessed through the zone mechanical energy active on command and through zone-based FISH functions zone.work.plastic.vol, zone.work.plastic.shear, zone.work.plastic.total, zone.work.elastic.vol, zone.work.elastic.shear, and zone.work.elastic.total.