Configurational forces

Derivation of configurational forces

Configurational body forces \((\vec{g}_{\textrm{body}})_{i}=g_{\textrm{body},i}\) for a hyper-elastic material are defined as [1]:

(1)\[g_{\textrm{body},i} := - \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} - \left( F^{\top} \right)_{ij} \cdot B_{j}\]

Note

The large displacement framework is used in this formulation of configurational forces. This has to be considered in Abaqus simulations by setting the flag NLGEOM=ON in order to compute the Cauchy stresses correctly.

The partial derivative of the Helmholtz energy density \(\Psi\) with respect to the undeformed coordinates \((\vec{X})_{i} = X_{i}\) is evaluated under a constant deformation gradient \((\textbf{F})_{ij} = F_{ij}\). This partial derivative is often called “internal configurational body forces”. The next term \(\left( F^{\top} \right)_{ij} \cdot B_{j}\) contains body forces \((\vec{B})_{j} = B_{j}\) and is often called “external configurational body forces”. Note, that \(\vec{B}(\vec{X})\) is defined on the undeformed configuration. Equation (1) cannot be computed easily, because the partial derivative is unknown. An alternative computation approach is shown in the following. For this we start with the Helmholtz energy density for a hyper-elastic material

(2)\[\Psi = \Psi \left( \vec{X}, \textbf{F} \left( \vec{X} \right) \right) \; ,\]

that depends on the undeformed coordinates \(\vec{X}\) and the deformation gradient \(\textbf{F}\). In this case the Helmholtz energy density \(\Psi\) is just the strain energy density. Next, we derive the Helmholtz energy density by \(\vec{X}\) using the chain rule.

(3)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} = \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} + \left( \frac{\partial \Psi}{\partial F_{jk}} \right)_{\vec{X}} \cdot \frac{\textrm{d} F_{jk}}{\textrm{d} X_{i}}\]

Note that \(\textrm{d} \Psi / \textrm{d} X_{i}\) and \(\textrm{d} F_{jk} / \textrm{d} X_{i}\) can be computed easily.

Piola-Kirchhoff stress

Furthermore, from FEM simulation we get the first Piola-Kirchhoff stress \((\textbf{P})_{ij} = P_{ij}\). For an elastic material in the large deformation framework, the first Piola-Kirchhoff stress

(4)\[P_{ij} = \left( \frac{\partial \Psi}{\partial F_{jk}} \right)_{\vec{X}}\]

corresponds to the partial derivative of the Helmholtz energy density \(\Psi\) with respect to the deformation gradient \(\textbf{F}\) at a fixed position \(\vec{X}\). [2]

Inserting equation (4) into equation (3) yields

(5)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} = \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} + P_{jk} \cdot \frac{\textrm{d} F_{jk}}{\textrm{d} X_{i}} \; .\]

Index swap of deformation gradient

The deformation gradient

(6)\[F_{jk} := \frac{\textrm{d} x_{j}}{\textrm{d} X_{k}}\]

is defined as the derivative of the deformed coordinates \((\vec{x})_{i} = x_{i} = X_{i} + U_{i}\) with respect to the undeformed coordinates \(\vec{X}\). This allows to swap the indices in the following equation

(7)\[\frac{\textrm{d} F_{jk}}{\textrm{d} X_{i}} = \frac{\textrm{d} }{\textrm{d} X_{i}} \frac{\textrm{d} x_{j}}{\textrm{d} X_{k}} = \frac{\textrm{d}^{2} x_{j} }{\textrm{d} X_{i} \; \textrm{d} X_{k}} = \frac{\textrm{d}^{2} x_{j} }{\textrm{d} X_{k} \; \textrm{d} X_{i}} = \frac{\textrm{d} }{\textrm{d} X_{k}} \frac{\textrm{d} x_{j}}{\textrm{d} X_{i}} = \frac{\textrm{d} F_{ji}}{\textrm{d} X_{k}} \; .\]

Inserting this in equation (5) yields

(8)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} = \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} + P_{jk} \cdot \frac{\textrm{d} F_{ji}}{\textrm{d} X_{k}} \; .\]

Product rule for derivative

For scalar functions \(f\) and \(g\), the product rule is well known as \((fg)'=f'g+fg'\). For the product of two tensors \(\textbf{F}\) and \(\textbf{P}\) the product rule is [3]:

(9)\[\frac{ \textrm{d} \left( \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{ \textrm{d} X_{k} } = \left( F^{\top} \right)_{ij} \cdot \frac{\textrm{d} P_{jk} }{ \textrm{d} X_{k} } + P_{jk} \cdot \frac{\textrm{d} F_{ji} }{ \textrm{d} X_{k} }\]

Note that the last term exists in equation (8). Pull out this term:

(10)\[P_{jk} \cdot \frac{\textrm{d} F_{ji} }{ \textrm{d} X_{k} } = \frac{ \textrm{d} \left( \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{ \textrm{d} X_{k} } - \left( F^{\top} \right)_{ij} \cdot \frac{\textrm{d} P_{jk} }{ \textrm{d} X_{k} }\]

And insert it into equation (8).

(11)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} = \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} + \frac{ \textrm{d} \left( \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{ \textrm{d} X_{k} } - \left( F^{\top} \right)_{ij} \cdot \frac{\textrm{d} P_{jk} }{ \textrm{d} X_{k} }\]

Body forces

The sum of all forces is zero. This is called the equilibrium of forces. For the static case this equilibrium is [2]:

(12)\[B_{j} = - \frac{\textrm{d} P_{jk}}{\textrm{d} X_{k}}\]

The body forces \(\vec{B}\) correspond to the divergence of the first Piola-Kirchhoff stress \(\textbf{P}\). Inserting this into equation (11) yields

(13)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} = \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} + \frac{ \textrm{d} \left( \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{ \textrm{d} X_{k} } + \left( F^{\top} \right)_{ij} \cdot B_{j} \; .\]

Configurational forces

We can rearrange equation (13)

(14)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} - \frac{ \textrm{d} \left( \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{ \textrm{d} X_{k} } = \left( \frac{\partial \Psi}{\partial X_{i}} \right)_{\textbf{F}} + \left( F^{\top} \right)_{ij} \cdot B_{j}\]

and compare the right side to the definition of configurational body forces in equation (1).

(15)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} - \frac{ \textrm{d} \left( \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{ \textrm{d} X_{k} } = -g_{\textrm{body}, i}\]

This equation contains only total derivatives of the Helmholtz energy density \(\Psi\), the deformation gradient \(\textbf{F}\) and the first Piola-Kirchhoff stress \(\textbf{P}\). These quantities can be computed from FEM results.

Prettify

Equation (15) can already be used. However, it is common to rearrange the equation into a more beautiful form with only one total derivative. For this, the derivative of the Helmholtz energy density is written as

(16)\[\frac{\textrm{d} \Psi}{\textrm{d} X_{i}} = \frac{\textrm{d} \Psi}{\textrm{d} X_{k}} \cdot \frac{\textrm{d} X_{k}}{\textrm{d} X_{i}} = \frac{\textrm{d} \Psi}{\textrm{d} X_{k}} \cdot \delta_{ik} = \frac{\textrm{d} \Psi \cdot \delta_{ik}}{\textrm{d} X_{k}} \; .\]

Insert this into equation (15).

(17)\[\frac{ \textrm{d} \left( \Psi \cdot \delta_{ik} - \left( F^{\top} \right)_{ij} \cdot P_{jk} \right) }{\textrm{d} X_{k}} = -g_{\textrm{body}, i}\]

Furthermore, Eshelby introduces the configurational stress (or energy momentum) tensor \((\textbf{CS})_{ik} = CS_{ik}\) [4]:

(18)\[CS_{ik} = \Psi \cdot \delta_{ik} - \left( F^{\top} \right)_{ij} \cdot P_{jk}\]

Using the configurational stress tensor, equation (17) is simplified to:

(19)\[\frac{ \textrm{d} CS_{ik} }{\textrm{d} X_{k}} = -g_{\textrm{body}, i}\]

Nodal configurational forces

For the computation of configurational forces from FEM results, it is common to compute nodal configurational forces. Each node is associated with a node volume by integrating the shape function of the i-th node \(\int_{\mathcal{B}} H_{i}(\vec{X})\, \textrm{d}V\). The sum of the node volumes is exactly the volume \(V=\int_{\mathcal{B}}1\,\textrm{d}V\) of the whole body \(\mathcal{B}\), because the sum of the shape functions is \(\sum_{i}H_{i}(\vec{X})=1\) at every position \(\vec{X}\). The same idea as used for the volume integral is also used for the nodal configurational forces. They are just weighted by the shape functions \((\vec{H})_{i} = H_{i}\) and are computed for each node. This is all done by our software. The nodal configurational force for the i-th node is [1]:

(20)\[g_{\textrm{nodal},ij} = \int_{\mathcal{B}} g_{\textrm{body},i} \cdot H_{j} \,\textrm{d}V\]

This can be written as:

(21)\[g_{\textrm{nodal},ij} = \int_{\mathcal{B}} -\frac{ \textrm{d} CS_{ik} }{\textrm{d} X_{k}} \cdot H_{j} \,\textrm{d}V\]

Integration by parts using Greens identity leads to:

(22)\[g_{\textrm{nodal},ij} = -\int_{\mathcal{\partial B}} CS_{ik} \cdot N_{k} \cdot H_{j} \,\textrm{d}A + \int_{\mathcal{B}} CS_{ik} \cdot \frac{ \textrm{d} H_{j} }{\textrm{d} X_{k}} \,\textrm{d}V\]

With the normal vector \((\vec{N})_{k}=N_{k}\) that points normal to the boundary \(\partial \mathcal{B}\) to the outside of the body \(\mathcal{B}\). Since the shape functions \(\vec{H}(\vec{X})\) are continuous analytical functions defined on the whole body \(\mathcal{B}\), it is easier to derive them instead of the configurational stresses \(\mathbf{CS}\), that are computed only at the integration points. Furthermore, the surface integral is zero, if the shape functions are all zero at the boundary. In this case the equation can be simplified to:

(23)\[g_{\textrm{nodal},ij} = \int_{\mathcal{B}} CS_{ik} \cdot \frac{ \textrm{d} H_{j} }{\textrm{d} X_{k}} \,\textrm{d}V\]

This integral is evaluated using the Gaussian integration that is commonly used in FEM. Our software conforce computes this integral and neglects the surface integral. Note, that neglecting the surface integral is not valid, for nodes lying at the boundary.

Displacement-based formulation

The displacement-based formulation (dbf) of the configurational forces replaces the deformation gradient \(\textbf{F}\) by the gradient of the displacements \((\vec{U})_{i} = U_{i}\).

Consequently, the dbf configurational stress is defined as [3]:

(24)\[CS_{\textrm{dbf}, ik} = \Psi \cdot \delta_{ik} - \frac{ \textrm{d} U_{j} }{\textrm{d} X_{i}} \cdot P_{jk}\]

Note, that there is a dependency between the deformation gradient \(\textbf{F}\) and the displacement gradient:

(25)\[F_{ik} := \frac{\textrm{d} x_{i}}{\textrm{d} X_{k}} = \delta_{ik} + \frac{\textrm{d} U_{i}}{\textrm{d} X_{k}}\]

An argument for dbf is, that it is numerically more stable for small displacements. Let’s for example consider a displacement gradient of \(1 \cdot 10^{-7}\). The deformation gradient adds one to this small numbers \(1 + 1 \cdot 10^{-7} = 1.0000001\). The computer might round this number to 1. Consequently, you remain more significant digits by using the displacement gradient instead of the deformation gradient.

The configurational stress \(\textbf{CS}\) is the dbf configurational stress \(\textbf{CS}_{\textrm{dbf}}\) minus the first Piola-Kirchhoff stress \(\mathbf{P}\).

(26)\[\begin{split}\begin{eqnarray} CS_{ik} & = & \Psi \cdot \delta_{ik} - \left( \delta_{ij} + \frac{\textrm{d} U_{j}}{\textrm{d} X_{i}} \right) \cdot P_{jk} \\ & = & \Psi \cdot \delta_{ik} - P_{ik} - \frac{\textrm{d} U_{j}}{\textrm{d} X_{i}} \cdot P_{jk} \\ & = & CS_{\textrm{dbf}, ik} - P_{ik} \end{eqnarray}\end{split}\]

In the absence of body forces \(\vec{B}\), the configurational body force \(\vec{g}_{\textrm{body}}\) corresponds to the dbf configurational body force \(\vec{g}_{\textrm{body, dbf}}\).

(27)\[\begin{split}\begin{eqnarray} g_{\textrm{body}, i} & = & - \frac{\textrm{d} CS_{ik}}{\textrm{d} X_{k}} \\ & = & - \frac{\textrm{d} (CS_{\textrm{dbf}, ik} - P_{ik})}{\textrm{d} X_{k}} \\ & = & - \frac{\textrm{d} CS_{\textrm{dbf}, ik} }{\textrm{d} X_{k}} + \frac{\textrm{d} P_{ik} }{\textrm{d} X_{k}} \\ & = & - \frac{\textrm{d} CS_{\textrm{dbf}, ik} }{\textrm{d} X_{k}} - B_{i} \\ & = & g_{\textrm{body,dbf}, i} - B_{i} \end{eqnarray}\end{split}\]

The nodal configurational forces can be converted using the following equation.

(28)\[\begin{split}\begin{eqnarray} g_{\textrm{nodal}, ij} & = & \int_{\mathcal{B}} g_{\textrm{body}, i} \cdot H_{j} \,\textrm{d}V \\ & = & \int_{\mathcal{B}} g_{\textrm{body, dbf}, i} \cdot H_{j} \,\textrm{d}V - \int_{\mathcal{B}} B_{i} \cdot H_{j} \,\textrm{d}V \\ & = & g_{\textrm{nodal,dbf}, ij} - B_{\textrm{nodal},ij} \end{eqnarray}\end{split}\]

Configurational forces in facture mechanics

Configurational forces are used to estimate the energy gradient with respect to a change in geometry. We call this gradient the energy release rate \(G=\partial \Pi / \partial \vec{p}\). The geometry change might be a movement of the crack tip position \(\vec{p}\). However, the definition of configurational forces in equation (1) does not derive by the crack tip position \(\vec{p}\) but derives by a coordinate \(\vec{X}\). This is not the same, as the following example illustrates. Figure 1 provides a graphical explanation.

Imagine you stand at a position \(\vec{X}\) and look at a point \(\vec{p}\). You can either move this point by \(\partial \vec{p}\) or you can move yourself by \(-\partial \vec{X}\). From your point of view, the point would look the same in both cases and you might guess that \(\partial \vec{p} = -\partial \vec{X}\). However, you also see the surrounding of the point. There are two cases:

1. In the first case, a second point \(\partial \vec{p'}\) lies in the surrounding of point \(\partial \vec{p}\). Figure 1 depicts the second point as red circle. If you move yourself by \(-\partial \vec{X}\), point \(\vec{p'}\) moves in common with point \(\vec{p}\) from your point of view. If you stand still and move point \(\vec{p}\) by \(\partial \vec{p}\) instead, the other point \(\vec{p'}\) stays fixed. Consequently, you can distinguish between a change of the point position and a change of the coordinates \(\partial \vec{p} \neq -\partial \vec{X}\).

2. In the second case, there is no other significant point in the surrounding of point \(\vec{p'}\). Consequently, \(\partial \vec{p} = -\partial \vec{X}\) holds true. We call this similarity.

Figure 1: Similarity of \(\textrm{d}X\) and \(-\textrm{d}p\)

Configurational forces are a valid estimation of the energy release rate \(G\) only in the second case. The following section explains how \(G\) is evaluated from the strain energy density \(\Psi\) using configurational forces. For a hyper-elastic material, the strain energy density \(\Psi\) depends on coordinates \(\vec{X}\), the deformation tensor \(\textbf{F}(\vec{X})\) and geometrical measures like the position of the crack tip \(\vec{p}\). In the absence of inertia and body forces like gravity \(\vec{B}=\vec{0}\), configurational body forces are gradients of \(\Psi\) with respect to coordinates \(\vec{X}\), but with a fixed deformation gradient \(\textbf{F}\) and a fixed position of the crack tip \(\vec{p}\).

(29)\[\vec{g}_{\textrm{body}}\left( \vec{X}, \textbf{F}, \vec{p} \right) = - \left( \frac{ \partial \Psi\left( \vec{X}, \textbf{F}, \vec{p} \right) }{ \partial \vec{X} } \right)_{ \textbf{F}, \vec{p} }\]

According to a similarity principle, the derivative by \(\partial \vec{X}\) can be replaced by the negative derivative by the position of the crack tip \(-\partial \vec{p}\).

(30)\[\vec{g}_{\textrm{body}}\left( \vec{X}, \textbf{F}, \vec{p} \right) = \left( \frac{ \partial \Psi\left( \vec{X}, \textbf{F}, \vec{p} \right) }{ \partial \vec{p} } \right)_{ \vec{X}, \textbf{F}}\]

In fracture mechanics, the strain energy \(\Pi\) is investigated instead of the energy density \(\Psi\). The strain energy is the integral over the energy density \(\Pi = \int \Psi \,\textrm{d}V\). Consequently, integrating the configurational body force over a certain crack-dominated region \(\mathcal{B}\) corresponds to the derivative of the strain energy.

(31)\[\left( \frac{ \partial \Pi \left( \mathcal{B}, \mathbf{F}, \vec{p} \right) }{ \partial \vec{p} } \right)_{\mathcal{B}, \mathbf{F}} = \int_{\mathcal{B}} \vec{g}_{\textrm{body}}\left( \vec{X}, \textbf{F}, \vec{p} \right) \,\textrm{d}V\]

Interpret this as a shift of a whole region \(\mathcal{B}\) instead of a single point. Note, that the region \(\mathcal{B}\) is fixed and configurational forces have to be zero at the boundary \(\partial \mathcal{B}\). The gradient can be decomposed into an energy release rate \(G\) and a unit vector \(\vec{v}\) with length one. The unit vector corresponds to the direction in which the crack will grow. This is the direction of the maximum energy dissipation. The energy release rate \(G\) can be compared to the fracture energy \(G_{c}\) to state whether the crack will grow or not.

(32)\[G \cdot \vec{v} = \int_{\mathcal{B}} \vec{g}_{\textrm{body}}\left( \vec{X}, \textbf{F}, \vec{p} \right) \,\textrm{d}V\]

In FEM, the evaluation of the integral can be simplified by using nodal configurational forces \(\vec{g}_{\textrm{nodal}, i}\) instead of the configurational body forces. In order to compute the energy release rate \(G\), all the user has to do is to call our function and then sum up the nodal configurational forces in the crack-dominated region and decompose the resulting nodal configurational force into the energy release rate \(G\) and the crack growth unit vector \(\vec{v}\).

(33)\[G \cdot \vec{v} = \sum_{i \in \mathcal{B}} \vec{g}_{\textrm{nodal}, i}\left( \textbf{F}, \vec{p} \right)\]

This demonstrates the advantage of configurational forces compared to the commonly used J-integral. Configurational forces provide the energy release rate as well as the crack growth direction, whereas the J-integral [5] only provides the energy release rate.

Plasticity

The implemented formulation does not support plasticity in general, since the Helmholtz energy density not only depends on \(X\) and \(F\), but also on plastic hardening parameters which are not considered here.

However, under the assumption of small strain plasticity, the formulation has already been used by Kolednik [6]. Kolednik proposes two modifications for the configurational stress:

• incremental plasticity, which considers only the elastic strain energy density (SENER)

\[CS_{\textrm{ep}, ik} = \Psi_{\textrm{elastic}} \cdot \delta_{ik} - (T^{\top})_{kj} \cdot P_{ij}\]

• deformation plasticity that considers both elastic and plastic strain energy densities (SENER+PENER)

\[CS_{\textrm{nlel}, ik} = (\Psi_{\textrm{elastic}} + \Psi_{\textrm{plastic}}) \cdot \delta_{ik} - (T^{\top})_{kj} \cdot P_{ij}\]

References

[1] (1,2)

R. Mueller and G. A. Maugin, “On material forces and finite element discretizations,” Computational Mechanics, vol. 29, no. 1, pp. 52–60, Jul. 2002, doi: 10.1007/s00466-002-0322-2.

[2] (1,2)

J. S. Bergstrom, “Mechanics of Solid Polymers: Theory and Computational Modeling”. Elsevier, 2015.

[3] (1,2)

M. E. Gurtin, Configurational forces as basic concepts of continuum physics. in Applied mathematical sciences, no. 137. New York: Springer, 2000.

[4]

J. D. Eshelby, “The force on an elastic singularity,” Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, vol. 244, no. 877, pp. 87–112, 1951.

[5]

J. R. Rice, “A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks,” Journal of Applied Mechanics, vol. 35, no. 2, pp. 379–386, Jun. 1968, doi: 10.1115/1.3601206.

[6]

O. Kolednik, R. Schöngrundner, and F. D. Fischer, “A new view on J-integrals in elastic–plastic materials,” Int J Fract, vol. 187, no. 1, pp. 77–107, May 2014, doi: 10.1007/s10704-013-9920-6.