15.2. Usage Guidelines

The following block command defines a CRK SECTION.

# Within Sierra scope...
BEGIN CRK SECTION <string>crk_name
  # Commands controlling approximation properties
  SHAPE FUNCTION TYPE = SPHERICAL|SPHERICAL_SURFACE_KRONECKER|
                        TWO_STAR|LINEAR_TET (SPHERICAL)
  SUPPORT SIZE = <real>support_size (1.1)
  SUPPORT SIZE TYPE = RELATIVE_LOCAL|RELATIVE_GLOBAL|
                      PHYSICAL (RELATIVE_LOCAL)
  BLOCK CONFORMING = ON|OFF (ON)

  # Commands controlling numerical integration
  FORMULATION = ELEMENT|NODAL_FROM_THEX|TET_FROM_THEX (ELEMENT)
  NUMBER OF INTEGRATION SUBCELLS = <int>num_sub_cells (1)

  # Controls which stress-strain conjugate pairs to use
  # ADAGIO uses a Cauchy-stress/rate-of-deformation conjugate
  # LOGARITHMIC uses Hencky strain/"Hencky" stress
  STRAIN MEASURE = ADAGIO|LOGARITHMIC (ADAGIO)

  # Controls when post processing occurs for performance
  POST PROCESS EVERY STEP = FALSE|TRUE (TRUE)

  # Creates face entities to apply contact to.  Required when
  # using Dash contact with CRK.
  CREATE CONTACT FACES = ON|OFF (OFF)

  # Commands controlling bond-based failure modeling in CRK
  FAILURE METHOD = BOND|NONE (NONE)
  FAILURE INDEX = <integer>failure_index
  FAILURE CRITERION = <real>bond_failure_criterion
  BOND OPENING DISPLACEMENT = <real>bond_opening_displacement
  BOND SECONDARY CRITERION = <real>bond_secondary_criterion
  SHOW BONDS = ON|OFF (OFF)

  # Controls which nodes are DOFs in the model
  # (for nearly-incompressible media)
  THEX ACTIVE NODES = NODES|EDGES|FACES|CENTROIDS (NODES)

END [CRK SECTION <string>crk_name]

The SHAPE FUNCTION TYPE, SUPPORT SIZE SUPPORT SIZE TYPE, and BLOCK CONFORMING options control the approximation properties of the CRK formulation.

The command SHAPE_FUNCTION_TYPE determines the kernel shape used in the CRK analysis. The option SPHERICAL is for the spherical-star where the elements that are partially or fully contained within a ball are included in the support. Similarly, the SPHERICAL_SURFACE_KRONECKER uses spherical-star kernels for the interior of the domain but with the support restricted to a one-star at surfaces, giving the same shape functions as a linear tetrahedral finite element on the surface and thus possessing the Kronecker delta property at these locations. The option TWO_STAR uses kernels with a two-star kernel definition. Finally, the option LINEAR_TET, gives linear tetrahedral finite element shape functions, equivalent to using a one-star kernel on a tetrahedral mesh.

Warning

The option SHAPE FUNCTION TYPE = TWO_STAR is only compatible with FORMULATION = ELEMENT described below. Additionally, the input mesh must be a non thex, tetrahedral topology.

The commands SUPPORT SIZE and SUPPORT SIZE TYPE control the support size and are only valid for SHAPE FUNCTION TYPE = SPHERICAL or SPHERICAL_SURFACE_KRONECKER. SUPPORT SIZE TYPE dictates whether the size is relative to the maximum length of attached edges (RELATIVE_LOCAL), relative to the global max element edge length (RELATIVE_GLOBAL), or absolute (PHYSICAL). The command SUPPORT SIZE sets the absolute or relative support size. Relative support sizes near but greater than 1 are recommended.

The command BLOCK CONFORMING = ON|OFF provides the option to make the kernel functions to conform to block boundaries or not. Non-conforming will give a behavior similar to traditional meshfree methods.

The command FORMULATION determines which type of numerical integration should be performed and supports three options. The option ELEMENT integrates using the elements of the input mesh as the strain-smoothing domains. NODAL_FROM_THEX specifies that nodal integration/strain-smoothing domains should be used and requires tetrahedral elements that are subdivided beforehand into hexahedral elements using the thex command in Cubit. The nodal integration domains are formed by aggregating all of the sub-tet hex elements connected to each node to form the barycentric dual, referred to as an integration “cell.” This results in the SCNI strategy [[1]]. The option TET_FROM_THEX integrates over the original tetrahedra before subdivision. For the latter two options, the sub-tet hex elements can be used for sub-cells in the \(\mathbf{\overline{F}}\) integration technique as described in Section 15.2.2.

If NODAL_FROM_THEX or TET_FROM_THEX are specified above, the nodes in the CRK domain will be divided into an “active” group and an “inactive” group. The “active” group will be the DOF-carrying nodes and usually represents the original tet nodes before the thex command is issued in Cubit. The “inactive” group will contain all other nodes which are effectively evaluation points used to setup the nodal integration domains used in CRK. These complementary sets may be adjusted using the THEX ACTIVE NODES command as described in Section 15.2.2.

Given its nonlocal nature, the CRK method typically lacks the Kronecker delta property and thus incurs additional post-processing steps for output purposes. The POST PROCESS EVERY STEP command controls whether post-processing procedures occur every step (TRUE) or during output steps only (FALSE). This improves overall performance as these computations are not necessary for the internal force routine and will decrease the total runtime if set to FALSE (contingent on the output frequency). However, derived output or user routines which depend on CRK variables on every step will not function correctly as a result. As such, the command POST PROCESS EVERY STEP = FALSE should only be used when user output, user functions, and other derived variables do not depend on CRK fields. The default is TRUE.

The command STRAIN MEASURE = ADAGIO|LOGARITHMIC controls which strain to use when evaluating the constitutive law. The default is ADAGIO which uses the rate of deformation (the default in Thecode). The option LOGARITHMIC uses the logarithmic strain, or Hencky strain, and its corresponding conjugate stress instead.

15.2.1. Contact

CRK currently supports Dash contact but the setup will be slightly different depending on whether the command FORMULATION is set to ELEMENT versus NODAL_FROM_THEX or TET_FROM_THEX. For both cases, the shape function type must be either SPHERICAL_SURFACE_KRONECKER or LINEAR_TET.

In the case that FORMULATION = ELEMENT, every node in the mesh will be labeled as an “active” node and Dash contact may be setup as if the simulation was using normal elements without any additional modification as given in the Sierra/SM User Manual.

In the case of FORMULATION = NODAL_FROM_THEX or TET_FROM_THEX, the CRK surface must be manually created using the command CREATE CONTACT FACES = ON. This command will skin the CRK block(s) and create a new block with the name <block_name>_contact_surface where <block_name> is the original block name. This new surface is necessary to setup CRK shape function weights and area densities that Dash requires. These weights spread the contact forces from inactive nodes to active ones associated with the problem unknowns. Once this surface is created, it may be referenced within the contact definition as, for example

# Within Sierra scope...
BEGIN CRK SECTION crk_section
  SHAPE FUNCTION TYPE = SPHERICAL_SURFACE_KRONECKER
  FORMULATION = NODAL_FROM_THEX
  CREATE CONTACT FACES = ON
END CRK SECTION crk_section
...
# Within Region scope...
BEGIN CONTACT DEFINITION collide
  SEARCH = Dash
  CONTACT SURFACE crk_surface CONTAINS block_1_contact_surface
  CONTACT SURFACE block_2 CONTAINS block_2

  BEGIN INTERACTION DEFAULTS
    GENERAL CONTACT = ON
    FRICTION MODEL = FRICTIONLESS
    INTERACTION BEHAVIOR = SLIDING
  END INTERACTION DEFAULTS
END CONTACT DEFINITION collide

where in the above block_1 is the CRK block and block_2 is an FEM block. Dash will then take the information present on the CRK surface including shape function values and weights in order to create the contact representation.

15.2.2. Advanced integration methods

By default, the \(\mathbf{\overline{F}}\) integration technique [[2]] is activated if the command FORMULATION = NODAL_FROM_THEX or TET_FROM_THEX and if the NUMBER OF INTEGRATION SUBCELLS is greater than one. In this case, the sub-tet hex elements will be used for sub-cells in the \(\mathbf{\overline{F}}\) method for improving performance with nearly incompressible materials. In this context a “sub-cell” is a collection of sub-tet hex elements which share an original tet node and represents a deviatoric domain of the integration cell. The number of sub-cells is specified with NUMBER OF INTEGRATION CELLS. Elements will be grouped into sub-cells using a k-means clustering algorithm if NUMBER OF INTEGRATION CELLS is less than the number of sub-tet hex elements for an integration cell. Each sub-tet hex element will be used as a sub-cell if NUMBER OF INTEGRATION CELLS is equal to or greater than the available elements for an integration cell. The volumetric portion of the deformation gradient is calculated over the cell. These are the nodal or tetrahedral domains for the NODAL_FROM_THEX and TET_FROM_THEX options, respectively.

The THEX ACTIVE NODES command controls which nodes should be active and defaults to NODES. This command is only valid with FORMULATION = NODAL_FROM_THEX or TET_FROM_THEX where the input mesh has been subdivided using the thex command in cubit. The thex command splits every tetrahedral element into four hexahedral elements, requiring the creation of nodes at the edge and face midpoints and the centroid of the tetrahedron. Any of the nodes in the resulting mesh can be activated (have degrees-of-freedom associated with them). The option NODES activates the nodes that were in the original tet mesh. The options EDGES, FACES and CENTROIDS activates the nodes that the thex process created on the input tetrahedral edges, faces and centroids, respectively. This option is useful for specifying additional “bubble” DOFs when using advanced integration methods such as \(\mathbf{\overline{F}}\) [[2]].

15.2.3. Bond-based fracture

Warning

The bond-based fracture method is only compatible with the command FORMULATION = NODAL_FROM_THEX.

CRK currently supports an experimental fracture formulation known as bond-based fracture. This formulation creates a “bond” between each node-node pair and maintains a damage variable associated with the bond which degrades the incremental deformation gradient and resulting internal force [[3]]. The failure criterion is specified using the FAILURE INDEX command to specify the emph{0-based} state variable index of the material model in use along with FAILURE CRITERION command which specifies the threshold the state variable must exceed before the node is considered “damaged.” When both nodes of a bond have reached the criterion, the bond is marked as damaged and begins to degrade linearly over a distance given by the BOND OPENING DISPLACEMENT command. This bond opening displacement is given by the relative displacement between the two nodes (in an abuse of terminology, the “stretch” of the bond).

The criterion above requires both nodes of a bond to be flagged as “damaged” before the bond can begin degrading. An alternative criterion may be specified via the command BOND SECONDARY CRITERION = <real>\(u_\textrm{second}\). For any bond that fails due to the first criterion, this option creates a failure plane passing through that bond’s midpoint normal to the bond. Any bond which shares a node with the first failed bond and intersects the failure plane may also fail due to a second criterion given by

\[u_{ab} - \bar{u}_{ab} \ge u_\textrm{second},\]

where \(u_{ab} = \lVert\mathbf{u}_b - \mathbf{u}_a\rVert\) is the relative displacement between nodes \(a\) and \(b\), and \(\bar{u}_{ab}\) is the relative displacement at the time one of the bond’s nodes becomes damaged due to the criterion given by FAILURE CRITERION and FAILURE INDEX described above. The bond will then begin to degrade over the BOND OPENING DISPLACEMENT as before.

Warning

The bonds may be viewed using the SHOW BONDS = ON command. However, this currently only works on a single processor.

15.2.4. Output

Note, CRK shape functions in general do not possess the Kronecker delta property, i.e.

\[N_I(x_J) \neq \delta_{IJ}.\]

As such, nodal variables such as displacement and velocity are not true displacements or velocities but rather the emph{generalized coefficients} of those fields. In addition, other quantities of interest will also be generalized (e.g. mass, force, etc.). These generalized values must be converted to their physical counterparts by interpolating these coefficients with their respective shape functions. These quantities are stored in the variables physical_displacement and physical_velocity. As a workaround, one can provide an alternative name in the results output to switch the meaning of displacement and physical_displacement as follows:

# Within Region scope...
BEGIN RESULTS OUTPUT my_output
  ...
  NODAL VARIABLES = physical_displacement as displacement
  NODAL VARIABLES = displacement as gen_displacement
  ...
END RESULTS OUTPUT my_output

which will output the displacement in a form amenable to Paraview or Ensight.

Warning

If an analysis contains both FEM and CRK formulations, then the displacement and velocity fields will be generalized for the CRK blocks and physical for the FEM blocks. Since physical_displacement does not exist in the FEM portions, the workaround mentioned above will output zero physical displacements for those blocks.

The output variables available for CRK are given below in Table 15.1 and Table 15.2 below. Bond-based fracture specific variables are given in Table 15.3.

Table 15.1 Nodal Variables for the CRK Formulation

Name	Type	Comments
physical_displacement	Vector_3D	The true displacement predicted by CRK
physical_velocity	Vector_3D	The true velocity predicted by CRK
displacement	Vector_3D	The generalized displacement predicted by CRK, i.e. the displacement coefficients
velocity	Vector_3D	The generalized velocity predicted by CRK, i.e. the velocity coefficients
nodal_status	Integer	Indicates whether a node is active (carrying a DOF) or inactive
max_edge_length	Real	Maximum edge length connected to the node of the original, non thex mesh
local_support_size	Real	Local support size radius when SHAPE FUNCTION TYPE is one of SPHERICAL or SPHERICAL_SURFACE_KRONECKER
num_neighbors	Integer	Number of active neighbors which have shape function value at the node
crk_contact_mass	Real	Generalized mass scattered to inactive nodes for Dash purposes

Table 15.2 Element Variables for the CRK Formulation

Name	Type	Comments
cell_id	Integer	The ID of the CRK strain smoothing cell
sub_cell_id	Integer	The ID of the CRK strain smoothing sub-cell
cell_volume	Real	The volume of the smoothing cell
sub_cell_volume	Real	The volume of the smoothing sub-cell

Table 15.3 Output Variables for Bond-Based FeFpmodelnameracture

Name	Type	Comments
node_damage	Nodal	Integer	Flag indicating if the failure criterion is met at a given node’s material point
damage_fraction	Element	Real	Fraction of damaged bonds over total bonds associated with an integration cell
bond_degradation	Element	Real	Scalar between 0 and 1 to indicate level of damage in the bond where 1 indicates no damage and 0 indicates fully damaged (only available when SHOW BONDS = ON)

---

[1]

Jiun-Shyan Chen, Cheng-Tang Wu, Sangpil Yoon, and Yang You. A stabilized conforming nodal integration for galerkin mesh-free methods. International Journal for Numerical Methods in Engineering, 50(2):435–466, 2001.

[2] (1,2)

Georgios Moutsanidis, Jacob Koester, Michael Tupek, Yuri Bazilevs, and Jiun-Shyan Chen. Treatment of near-incompressibility in meshfree and immersed-particle methods. Computational Particle Mechanics, 2019.

[3]

Jacob Koester, Michael Tupek, and Scott Mitchell. An agile design-to-simulation workflow using a new conforming moving least squares method. Technical Report SAND2019-11851, Sandia National Laboratories, Albuquerque, NM, 2019.