BilinearMixedModeCohesiveZoneModel

Description

The BilinearMixedModeCohesiveZoneModel object computes a global traction vector from a local traction vector that is calculated from a bilinear mixed mode cohesive zone model. Theoretical equations and their application to particle fuel is described and discussed in Jiang et al. (2021).

This object inherits from PenaltySimpleCohesiveZoneModel and can include penalty mortar mechanical contact.

This object implements the bilinear mixed mode traction separation law using the nodal-based mortar quantities that are employed in the MOOSE formulation. Such quantities are identified with a tilde at a secondary surface node . The definitions of mortar surfaces and mortar quantities are identical to those presented in previous technical reports (Martin Recuero and Yushu (2022) and Recuero et al. (2022)), which are in turn based on existing discretization approaches (Wohlmuth (2011) and Puso et al. (2008)).

Quadratic failure criterion

The softening of the interface in the bilinear mixed mode traction separation law is given by a quadratic failure criterion. The onset of softening of the interface takes place when the traction state reaches the following equality: (1) where , , and are the (positive) normal traction, the shear traction along direction one, and the shear traction along direction two, respectively. These three traction components are computed at each node using interpolated, weighted quantities. , , and are the interlaminar tensile and shear strengths, which are model parameters.

Mode mixity

The displacement on the interface can be split into normal and tangential parts. Therefore, the mixed mode relative displacement can be expressed as: (2) where subscripts , and and denote the normal direction and the two tangential directions, respectively.

The formulation assumes that the interface traction generated by the displacements is proportional; in essence, a traction component is equal to a penalty factor times the displacement jump component (, where is the penalty factor). If , the mixity ratio is obtained as: (3)

Assuming , the uni-mode relative displacements at the beginning of softening (denoted by superscript 0'') are given by: (4)

The mixed mode relative displacement corresponding to the initiation of softening, , is defined by the following multi-part definition: (5)

Mode mixity Delamination propagation prediction

The new CZM mortar implementation also accounts for two types of delamination propagation models. (1) A power law-based criterion and (2) The Benzeggagh-Kenane (BK) criterion (see, e.g.,Benzeggagh and Kenane (1996) and Bui (2011)).

Power law-based delamination criterion

The power law criterion may be described as: (6) where and are the critical energy release rates for mixed modes; whereas and are the critical energy release rates for pure modes. The power exponent is can be obtained experimentally.

During the delamination or debonding process, a relative interface displacement may be reached such that cohesive fores completely disappear, i.e. the cohesive model ceases to influence boundary tractions. Such a mixed-mode displacement value is defined as follows (note that we carry the node-wise definition, ): (7)

Benzeggagh-Kenane (BK) delamination criterion

The BK criterion can be written as (Benzeggagh and Kenane (1996)): (8) where , is the total critical strain energy release rate, and is a model parameter.

The mixed-mode total decohesion displacements are, for the BK criterion: (9)

Constitutive model

The CZM constitutive models relate the weighted interface jump components to the interface traction in the local frame of reference. The traction, in the local frame of reference, can be expressed as: (10) where is the component of the local traction and is an effective stiffness. The constitutive components are computed as: (11) where represents the CZM weighted interface jump with shear components, i.e. it applies shear stiffness, and normal stiffness only if a gap exists. applies a normal traction that prevents interpenetration. A few notes about the previous equation:

- The action of can potentially be replaced with (mortar) mechanical contact enforcement.
- The local traction is computed nodally with weighted quantities, but is interpolated and applied on the primary and secondary surfaces using Lagrange interpolation functions.
- The constitutive model can be expanded as needed in the user object.

The damage, , is computed as follows (we omit the nodal subscript):

Viscous regularization

Sudden degradation of stiffness may be give rise to the appearance of large displacements. To alleviate this problematic numerical behavior in quasti-static simulations, we include in the constitutive equations a viscous damage instead:

Example of usage:

[UserObjects]
  [czm_uo]
    type = BilinearMixedModeCohesiveZoneModel
    primary_boundary = 'top_bottom'
    secondary_boundary = 'bottom_top'
    primary_subdomain = 10000
    secondary_subdomain = 10001

    disp_x = disp_x
    disp_y = disp_y
    friction_coefficient = 0.1 # with 2.0 works
    secondary_variable = disp_x

    penalty = 0e6
    penalty_friction = 1e4
    use_physical_gap = true

    correct_edge_dropping = true
    normal_strength = N
    shear_strength = 1e3
    viscosity = 1e-3
    penalty_stiffness = 1e6

    power_law_parameter = 2.2
    GI_c = 1e3
    GII_c = 1e2
    displacements = 'disp_x disp_y'
  []
[]
(modules/contact/test/tests/cohesive_zone_model/mortar_czm.i)

Input Parameters

  • disp_xThe x displacement variable

    C++ Type:std::vector<VariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The x displacement variable

  • disp_yThe y displacement variable

    C++ Type:std::vector<VariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The y displacement variable

  • displacementsThe string of displacements suitable for the problem statement

    C++ Type:std::vector<VariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The string of displacements suitable for the problem statement

  • friction_coefficientThe friction coefficient ruling Coulomb friction equations.

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The friction coefficient ruling Coulomb friction equations.

  • penaltyThe penalty factor

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The penalty factor

  • primary_boundaryThe name of the primary boundary sideset.

    C++ Type:BoundaryName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the primary boundary sideset.

  • primary_subdomainThe name of the primary subdomain.

    C++ Type:SubdomainName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the primary subdomain.

  • secondary_boundaryThe name of the secondary boundary sideset.

    C++ Type:BoundaryName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the secondary boundary sideset.

  • secondary_subdomainThe name of the secondary subdomain.

    C++ Type:SubdomainName

    Unit:(no unit assumed)

    Controllable:No

    Description:The name of the secondary subdomain.

  • secondary_variablePrimal variable on secondary surface.

    C++ Type:VariableName

    Unit:(no unit assumed)

    Controllable:No

    Description:Primal variable on secondary surface.

Required Parameters

  • adaptivity_penalty_normalSIMPLEThe augmented Lagrange update strategy used on the normal penalty coefficient.

    Default:SIMPLE

    C++ Type:MooseEnum

    Unit:(no unit assumed)

    Options:SIMPLE, BUSSETTA

    Controllable:No

    Description:The augmented Lagrange update strategy used on the normal penalty coefficient.

  • aux_lmAuxiliary variable that is utilized together with the penalty approach to interpolate the resulting contact tractions using dual bases.

    C++ Type:std::vector<VariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:Auxiliary variable that is utilized together with the penalty approach to interpolate the resulting contact tractions using dual bases.

  • correct_edge_droppingFalseWhether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether to enable correct edge dropping treatment for mortar constraints. When disabled any Lagrange Multiplier degree of freedom on a secondary element without full primary contributions will be set (strongly) to 0.

  • debug_meshFalseWhether this constraint is going to enable mortar segment mesh debug information. An exodusfile will be generated if the user sets this flag to true

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether this constraint is going to enable mortar segment mesh debug information. An exodusfile will be generated if the user sets this flag to true

  • disp_zThe z displacement variable

    C++ Type:std::vector<VariableName>

    Unit:(no unit assumed)

    Controllable:No

    Description:The z displacement variable

  • ghost_higher_d_neighborsFalseWhether we should ghost higher-dimensional neighbors. This is necessary when we are doing second order mortar with finite volume primal variables, because in order for the method to be second order we must use cell gradients, which couples in the neighbor cells.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether we should ghost higher-dimensional neighbors. This is necessary when we are doing second order mortar with finite volume primal variables, because in order for the method to be second order we must use cell gradients, which couples in the neighbor cells.

  • ghost_point_neighborsFalseWhether we should ghost point neighbors of secondary face elements, and consequently also their mortar interface couples.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether we should ghost point neighbors of secondary face elements, and consequently also their mortar interface couples.

  • interpolate_normalsFalseWhether to interpolate the nodal normals (e.g. classic idea of evaluating field at quadrature points). If this is set to false, then non-interpolated nodal normals will be used, and then the _normals member should be indexed with _i instead of _qp

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether to interpolate the nodal normals (e.g. classic idea of evaluating field at quadrature points). If this is set to false, then non-interpolated nodal normals will be used, and then the _normals member should be indexed with _i instead of _qp

  • minimum_projection_angle40Parameter to control which angle (in degrees) is admissible for the creation of mortar segments. If set to a value close to zero, very oblique projections are allowed, which can result in mortar segments solving physics not meaningfully, and overprojection of primary nodes onto the mortar segment mesh in extreme cases. This parameter is mostly intended for mortar mesh debugging purposes in two dimensions.

    Default:40

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:Parameter to control which angle (in degrees) is admissible for the creation of mortar segments. If set to a value close to zero, very oblique projections are allowed, which can result in mortar segments solving physics not meaningfully, and overprojection of primary nodes onto the mortar segment mesh in extreme cases. This parameter is mostly intended for mortar mesh debugging purposes in two dimensions.

  • penalty_frictionThe penalty factor for frictional interaction. If not provide, the normal penalty factor is also used for the frictional problem.

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The penalty factor for frictional interaction. If not provide, the normal penalty factor is also used for the frictional problem.

  • periodicFalseWhether this constraint is going to be used to enforce a periodic condition. This has the effect of changing the normals vector for projection from outward to inward facing

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether this constraint is going to be used to enforce a periodic condition. This has the effect of changing the normals vector for projection from outward to inward facing

  • primary_variablePrimal variable on primary surface. If this parameter is not provided then the primary variable will be initialized to the secondary variable

    C++ Type:VariableName

    Unit:(no unit assumed)

    Controllable:No

    Description:Primal variable on primary surface. If this parameter is not provided then the primary variable will be initialized to the secondary variable

  • prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.

    C++ Type:MaterialPropertyName

    Unit:(no unit assumed)

    Controllable:No

    Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.

  • use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

  • use_physical_gapFalseWhether to use the physical normal gap (not scaled by mortar integration) and normalize the penalty coefficient with a representative lower-dimensional volume assigned to the node in the contacting boundary. This parameter is defaulted to 'true' in the contact action.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether to use the physical normal gap (not scaled by mortar integration) and normalize the penalty coefficient with a representative lower-dimensional volume assigned to the node in the contacting boundary. This parameter is defaulted to 'true' in the contact action.

Optional Parameters

  • GII_cCritical energy release rate in shear direction.

    C++ Type:MaterialPropertyName

    Unit:(no unit assumed)

    Controllable:No

    Description:Critical energy release rate in shear direction.

  • GI_cCritical energy release rate in normal direction.

    C++ Type:MaterialPropertyName

    Unit:(no unit assumed)

    Controllable:No

    Description:Critical energy release rate in normal direction.

  • lag_displacement_jumpFalseWhether to use old displacement jumps to compute the effective displacement jump.

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether to use old displacement jumps to compute the effective displacement jump.

  • mixed_mode_criterionBKOption for mixed mode propagation criterion.

    Default:BK

    C++ Type:MooseEnum

    Unit:(no unit assumed)

    Options:POWER_LAW, BK

    Controllable:No

    Description:Option for mixed mode propagation criterion.

  • normal_strengthTensile strength in normal direction.

    C++ Type:MaterialPropertyName

    Unit:(no unit assumed)

    Controllable:No

    Description:Tensile strength in normal direction.

  • penalty_stiffnessPenalty stiffness for CZM.

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:Penalty stiffness for CZM.

  • power_law_parameterThe power law parameter.

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:The power law parameter.

  • regularization_alpha1e-10Regularization parameter for the Macaulay bracket.

    Default:1e-10

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:Regularization parameter for the Macaulay bracket.

  • shear_strengthTensile strength in shear direction.

    C++ Type:MaterialPropertyName

    Unit:(no unit assumed)

    Controllable:No

    Description:Tensile strength in shear direction.

  • viscosity0Viscosity for damage model.

    Default:0

    C++ Type:double

    Unit:(no unit assumed)

    Controllable:No

    Description:Viscosity for damage model.

Bilinear Mixed Mode Traction Parameters

  • allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).

  • control_tagsAdds user-defined labels for accessing object parameters via control logic.

    C++ Type:std::vector<std::string>

    Unit:(no unit assumed)

    Controllable:No

    Description:Adds user-defined labels for accessing object parameters via control logic.

  • enableTrueSet the enabled status of the MooseObject.

    Default:True

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:Yes

    Description:Set the enabled status of the MooseObject.

  • execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.

    Default:0

    C++ Type:int

    Unit:(no unit assumed)

    Controllable:No

    Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.

  • force_postauxFalseForces the UserObject to be executed in POSTAUX

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Forces the UserObject to be executed in POSTAUX

  • force_preauxFalseForces the UserObject to be executed in PREAUX

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Forces the UserObject to be executed in PREAUX

  • force_preicFalseForces the UserObject to be executed in PREIC during initial setup

    Default:False

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Forces the UserObject to be executed in PREIC during initial setup

  • use_displaced_meshTrueWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

    Default:True

    C++ Type:bool

    Unit:(no unit assumed)

    Controllable:No

    Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

Input Files

References

  1. ML Benzeggagh and MJCS Kenane. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Composites science and technology, 56(4):439–449, 1996.[BibTeX]
  2. Q.V. Bui. A modified benzeggagh-kenane fracture criterion for mixed-mode delamination. Journal of Composite Materials, 45(4):389–413, 2011.[BibTeX]
  3. Wen Jiang, Gyanender Singh, Jason D Hales, Aysenur Toptan, Benjamin W Spencer, and Stephen R Novascone. Efficient failure probability calculations and modeling interface debonding in TRISO particles with BISON. Technical Report, Idaho National Lab.(INL), Idaho Falls, ID (United States), 2021.[BibTeX]
  4. Antonio Martin Recuero and Dewen Yushu. Implement and test 3D mortar contact in bison. Technical Report, Idaho National Laboratory (INL), Idaho Falls, ID (United States), 2022.[BibTeX]
  5. Michael A Puso, TA Laursen, and Jerome Solberg. A segment-to-segment mortar contact method for quadratic elements and large deformations. Computer Methods in Applied Mechanics and Engineering, 197(6-8):555–566, 2008.[BibTeX]
  6. Antonio Recuero, Alexander Lindsay, Dewen Yushu, John W Peterson, and Benjamin Spencer. A mortar thermomechanical contact computational framework for nuclear fuel performance simulation. Nuclear Engineering and Design, 394:111808, 2022.[BibTeX]
  7. Barbara Wohlmuth. Variationally consistent discretization schemes and numerical algorithms for contact problems. Acta Numerica, 20:569–734, 2011.[BibTeX]