Neper2CAE and PyCiGen: Scripts to generate polycrystals and interface elements in Abaqus

Two codes have been developed: one allows columnar polycrystal geometries generated by software such as Neper or DREAM3D to be imported into the Complete Abaqus Environment for further geometric manipulation and meshing with arbitrary element types; the other code generates zero thickness, four noded interface elements at the boundaries between the crystals. Together these allow intergranular fracture simulations to be performed easily in Abaqus. This approach can automatically create different grain structures and allows the mechanical properties of metals and other crystalline materials to be simulated. This is important as it enables the sample to sample variability of fracture and mechanical properties that is observed in experiments to be reproduced. These scripts extend the simulation capabilities of the finite element solver Abaqus.


Motivation and significance
To simulate the mechanical properties of metals requires analysis of the elastic and plastic response of the microstructure [1]. Metals consist of grains with a crystalline structure, whose dimensions can vary from nanometres to millimetres. Each grain has a specific shape and orientation. The simulation of mechanical tests can be carried out using finite element solvers, such as Abaqus [2]. User material (UMAT) subroutines can be used to model the mechanical properties of different grains [3,4]. Cohesive interfaces implemented through a zero thickness user element (UEL) subroutine are used to model fracture nucleation and propagation [5][6][7].
Software for the generation of polycrystal geometries are available, such as Neper [8] and DREAM3D [9]. They can generate polycrystals with arbitrary grain sizes, distributions and orientations. They can also import experimental images. Meshing capabilities of the polycrystal geometry are usually included in the above-mentioned software. However, the geometric shapes that can be generated and meshing options are limited. Specifically, a graphical user interface that allows a wide variety of geometrical manipulation options is not available.
Usually, only parallelepiped geometries with four node tetrahedral elements can be generated. Hexahedral elements are preferred to tetrahedral ones, which give an excessively stiff response during bending [10]. Four node tetrahedral elements exhibit the so-called volumetric locking [11]. This numerical phenomenon takes place when the displacement field has nearly zero volumetric strain. The shape functions of tetrahedral elements are less accurate in capturing those deformations and this results in pairs of neighbouring elements carrying positive and negative volumetric strains respectively. Therefore, spurious high stresses are introduced. This is particularly true for crystal plasticity constitutive models, in which the plastic deformation is modelled as an isochoric process. [12][13][14].
Abaqus CAE (Complete Abaqus Environment) offers more advanced capabilities. The geometry can be modified using extrusions of arbitrary shapes and multiple parts can be assembled in a model, therefore samples with arbitrary shapes can be reproduced. Moreover, internal voids can be introduced. If external devices apply load to the sample, they can be included in the model. Arbitrary boundary conditions and contact between parts can be introduced. The mesh generation features of Abaqus CAE are also more advanced. Different element types and meshing algorithms can be used.
Therefore, linking polycrystal generation software, such as Neper, with Abaqus CAE represents a strong need for the scientific and engineering community. Such a link is provided by the program Neper2CAE, described in this manuscript.
Automated generation of polycrystal geometries is also useful when running a set of simulations with different grain size and shapes. Thus, the sample to sample variability of fracture and mechanical properties can be understood and safety limits can be established during component design.
Grain boundaries, which are the 2D interfaces between neighbouring grains, provide crack nucleation and propagation sites in a range of metals [15]. This phenomenon can be modelled using interface elements, whose mechanical properties are described by a traction-separation law [16]. Interface elements constitute a 2D mesh on the grain boundaries. However, commercial finite element software does not offer the possibility to automatically generate zero thickness cohesive elements. Available software for polycrystal generation can create only triangular interface elements. If hexahedral elements are used inside the grains, quadrilateral interface elements must be used. Therefore, an additional program is needed, which is PyCiGen, described in this manuscript.
The software has been developed specifically for polycrystalline metals, but its applicability can be easily extended to different materials with interfaces, e.g. composite materials.

Software description: polycrystal generation
In this section, the main Abaqus scripting commands used for the polycrystal generation are reported. A Neper parser (.tess file) is also present in the code developed, but it is not described in detail in the manuscript and the reader is referred to the code repository.
A 2D polycrystal is generated using Neper [8], as shown in Fig. 1. The generated .tess file contains the coordinates (x, y) of the vertices v i of the grains after the **vertex keyword. These are introduced in the upper surface of the Abaqus geometry using the command: where z is the user-defined depth of the parallelepiped geometry along the z axis.
The edges that connect two vertices v i and v j of the grains are described after the keyword **edge in the .tess file. They are included in the Abaqus geometry as partitions of the upper surface using the command:

faces=pickedFaces)
A 2D grain structure is therefore generated on the upper surface, as shown by the grain inside the red dashed square in Fig. 1.
The geometry is partitioned into columnar grains by sweeping the edges e i , generated on the upper surface, along a sweep direction, which correspond to the negative z direction, as shown in Fig. 1. This is done using the command: At this point, one cell for each grain is present. A user material is defined using: The mechanical constants indicate the orientation of the grain and the material type. In this way, dual-phase materials, with grains constituted of different crystal structures, can be included. The corresponding solid section is created and assigned to each grain using: Finally the Abaqus model is ready and can be transformed using the Abaqus CAE capabilities. Boundary conditions can be assigned and a mesh can be created. Usage instructions: in the file main.py, the following parameters must be set: width, height and depth. These will be the dimensions of the parallelepiped geometry along the x, y and z axes generated in Abaqus. The parameter determining the number of grains Ngrains must also be set. The code is then executed by: python3.6 main.py

Software description: interface elements generation
In this section, a mesh of hexahedral elements is considered where each bulk element has 8 nodes, corresponding to the vertices of the hexahedron. Adjacent elements are connected by a face, which contains 4 nodes. Each node and each element is characterised by a unique positive integer index.
First a data structure called ''connectivity'' is constructed that stores the indices of all bulk elements that contain a particular node n i . Bulk element sets are built by Abaqus that contain all the elements belonging to a certain grain.
Interface elements are built at the interface between grains. Therefore, the nodes that are at those interfaces have to be identified. For instance, in Fig. 2, the node n 1 belongs to three different grains (G1, G2 and G3), while node n 2 belongs to two different grains (G1 and G2).
In order to identify interface nodes, the ''connectivity'' is used: given a node n i , the bulk elements that contain that node and the grains to which they belong can be found. If the bulk elements that contain n i belong to different grains, then n i is an interface node.
Adjacent bulk elements at the grain boundary are identified as elements containing interface nodes. ''Face'' objects are created at the interface between two adjacent bulk elements. The pairs of adjacent bulk elements are identified by checking if they share four interface nodes. In Fig. 2, each ''face'' object is represented by four dashed lines connecting nodes of bulk elements in different grains.
Interface nodes are multiplied. For instance, in Fig. 2, the node n 1 is triplicated because it belongs to three different grains. Two new nodes n ′ 1 and n ′′ 1 are created. The node n 2 is duplicated and a new node n ′ 2 is created. Each duplicated node contains a variable that points to the original node. Each duplicated node is assigned to a different bulk element belonging to a different grain. For instance, in Fig. 2, n 1 , n ′ 1 and n ′′ 1 are assigned to three different grains. After this operation, there are no more interface nodes: each node belongs to a unique grain.
At this point, the face object contains information about the 8 nodes that will form the interface element and information about the 2 adjacent bulk elements. The first 4 nodes (i 1 , i 2 , i 3 , i 4 ) of the interface element belong to one grain and the others (i 5 , i 6 , i 7 , i 8 ) belong to the other grain.
In a finite element simulation, a deformation vector is associated to each node and a traction-separation law is used, which is based on that deformation. Therefore, the order of the nodes in the interface element is important. At zero deformation, node i 1 coincides with i 5 , node i 2 with i 6 , node i 3 with i 7 , node i 4 with i 8 , as shown in Fig. 3(b). The distances between these pairs of nodes determine the interface separation.
A right-handed reference system XYZ is built, as shown in Fig. 3(b). X points from i 1 to i 2 , Y is perpendicular to X and belongs to the plane containing the nodes (i 1 , i 2 , i 3 , i 4 ). The Z axis points from grain 2 (G2) to grain 1 (G1) in Fig. 3(b). The direction of Z is important for the finite element solver because it discriminates between opening and closing of the cohesive interface. In Fig. 3(b), the quadrilateral constituted by the nodes (i 5 , i 6 , i 7 , i 8 ) has moved towards the positive Z direction, indicating that the cohesive interface is opening. Z must always point from the grain that contains the interface nodes (i 1 , i 2 , i 3 , i 4 ) to the grain that contains the interface nodes (i 5 , i 6 , i 7 , i 8 ).
To determine the correct direction of Z and the correct reference system XYZ , it is important to associate properly the  nodes of the interface element (i 1 , i 2 , i 3 , i 4 ) with the nodes of the corresponding bulk element. The criterion is based on the righthand rule: if the right hand fingers curl along the path i 1 → i 2 → i 3 → i 4 , the thumb must point towards the other grain.
This direction determines the Z axis. For the other adjacent bulk element if the fingers curl along the path i 5 → i 6 → i 7 → i 8 , the right thumb must also point along Z towards the grain containing This association is done by using the following property of the hexahedral elements generated by Abaqus: the elements can be arbitrarily oriented in space but the node numbering is always the one represented in Fig. 3(a). Six possible faces could be associated with interface elements. The face object contains information about the adjacent bulk elements. Therefore, a ''facetype'' variable is defined that contains the information about the bulk element face that lies on the grain boundary.
For instance, if a bulk element face with nodes (n 1 , n 2 , n 3 , n 4 ) is associated with the interface element nodes (i 1 , i 2 , i 3 , i 4 ), the right-hand rule explained above applies. The correspondence between nodes of the bulk elements and nodes of the interface elements will be: i 1 ≡ n 4 , i 2 ≡ n 3 , i 3 ≡ n 2 , i 4 ≡ n 1 . If the same bulk element face is associated with the interface element nodes (i 5 , i 6 , i 7 , i 8 ), then the correspondence is: i 5 ≡ n 1 , i 6 ≡ n 2 , i 7 ≡ n 3 , i 8 ≡ n 4 . This is shown in Fig. 3(b) for the bulk element in grain 1 (G1).
For the six face-types, the correspondence between nodes of the interface elements and nodes of the bulk elements is reported in Table 1.
The behaviour of the interface elements has been studied by simulating two bulk elements and one interface element at the common interface between them. Fig. 4(a) shows the result of this test. The plot shows the traction on the interface element as Table 1 Possible correspondences between bulk element nodes and interface element nodes.

Conclusions
Two programs have been developed to generate polycrystal representative volumes and interface elements at the grain boundaries. These features are not available in Abaqus CAE but they are normally required to model mechanical properties at the micrometre length scale.
The algorithm for interface element generation is of general interest. It represents an example of interface detection in a mesh. The strategy used to find the orientation of the interface elements, described in Section 3, is based on the ordering of the interface element nodes and on the specific correspondence between interface element nodes and bulk element nodes. It can be used to find the directions pointing towards the inside and outside of a generic element set.
This software has been used in a recent study on the intergranular fracture of cast α-uranium [17,18] and it is available in the following repositories [19,20]. The implementation of different types of interface elements can be a future improvement of PyCiGen.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.