Determination of the effective 2D material stiffness matrix using ABAQUS

A cross‐scale homogenization approach for multi‐phase materials is to calculate an effective material stiffness in the smaller scale for the multiple states. These detailed results are used to identify a calculation rule for the effective material stiffness in the larger scale.

A cross-scale homogenization approach for multi-phase materials is to calculate an effective material stiffness in the smaller scale for the multiple states. These detailed results are used to identify a calculation rule for the effective material stiffness in the larger scale.
The calculation of the stiffness matrix using ABAQUS for 3D finite elements was already described by [1]. This can be transferred to the calculation of the 2D stiffness matrix in R 3×3 without much effort. However, some 2D elements in ABAQUS require the specification of a matrix in R 4×4 . In this article, the methods for calculating this type of stiffness matrix in 2D are presented. Goldberg et al. described in [1,2] how to determine the effective 3D material stiffness matrix in ABAQUS. The basis for this is a representative volume element (RVE) in which the displacements of the opposing surfaces are coupled with each other via periodic boundary conditions, whereby the Hill-Mandel conditions are fulfilled.
For the determination of the 2D material stiffness matrix in R 3×3 this procedure, including the equations, can be adapted to the 2D formulation. In Fig. 1 a 2D RVE with three pilot nodes is shown as an example. One pilot node (P 0 ) is used to suppress the rigid body movements, the other two pilot nodes (P 1 , P 2 ) are used to apply the displacement boundary conditions.
The reaction forces Q Pi are calculated by ABAQUS at the pilot nodes for each spatial direction. From these and the given displacements u Pi the effective Cauchy stress σ can be calculated (compare [1]): Here l i e i = OP i − OP 0 , A is the area in the undeformed state and I is the identity matrix.
In general, the materials show a non-linear behavior, i.e. a material stiffness that changes with deformation, which can be represented by nonlinear constitutive equations. To calculate the effective material stiffness, the linear HOOKE'S constitutive law σ = C ·· ε is used, which corresponds to a linearizaiton in the undeformed initial configuration.
The calculation of the stiffness matrix depends on three different DIRICHLET boundary conditions, as listed in Eq. 2. Here, every displacement of the pilot nodes ia applied by a sufficiently small scalar λ. The reaction forces at the pilot nodes calculated by ABAQUS are used to determine the stress by using Eq. 1. Afterwards, each entry of the stiffness matrix can be calculated directly from the stress response using Eq. 3. After the computation, the stiffness matrix should be symmetrized to minimize numerical errors.   2 Determination of the material stiffness matrix in R 4×4 When using plane strain in ABAQUS, the UMAT requires a material stiffness matrix in R 4×4 for the CPE4 element as default, which is used for stress calculation as shown in Eq. 4. To calculate this material stiffness matrix, the same ABAQUS model, as shown in Sec. 1, is used. The displacement boundary condition is again applied via the pilot nodes. However, the stress component σ 33 is now required for the calculation of the stiffness matrix. This stress cannot be calculated from the reaction forces at the pilot nodes, because in ABAQUS only the reaction forces in the in-plane directions are provided. Therefore, to determine the effective material stiffness, the stresses at the GAUSSIAN points are calculated and integrated to a total stress over the surface. The calculation rules for this are shown in Eq. 5, where the index G indicates values at the GAUSSIAN points.
Subsequently, the following equation results for the effective stiffness matrix under consideration of the symmetry: The value of C33 can be chosen arbitrarily. Here it was chosen to ensure that the matrix will not become ill-conditioned. As before, the matrix will be symmetrized.

Example
A square area with a fiber reinforced material is used as numerical example. Therefore, a 2D RVE is modeled with different fiber directions. The material model is based on a nonlinear NEO-HOOKE model with an additional fiber reinforced material part, compare Eq. 7 (see also [3]). The fiber direction is given by a.
The exact values of the material parameter G(j), K(j) and G F depend on whether the element belongs to the matrix or fiber. In Fig. 2 a configuration with horizontal fiber arrangement is shown and the determined corresponding material stiffness matrix. The stiffening effect of the fiber can be recognized by the significantly larger value of the component C 11 in the determined stiffness matrix. Similarly, when using an angle of 26.6 • to the horizontal as fiber direction (see Fig. 3) the stiffness is correctly determined with this method. In the horizontal direction, a slight stronger stiffening is visible in the C 11 component and concerning shear, a distinct stiffening is visible in the C 14 , C 24 , C 34 components.