Cohesive Zone Modeling

Cohesive zone modeling can be used to model adhesive and bonded interfaces and corresponding crack initiation and propagation. There are multiple methods using adhesive and bonded interfaces which can currently be modeled in OptiStruct.

Depending on the method, thickness of the adhesive/bonded interface can be defined (for damage model-based cohesive elements). The locations where cracks are expected should be identified and defined as cohesive zones.

Adhesive/Bonded Interface Implementation Methods

Cohesive zones are used to model Adhesive/Bonded Interfaces, which are typically locations where parts are connected with a glue or glue-like material.

There are currently two approaches for simulating such interfaces:

Traction-Opening Curves Method
Damage Models Method

Traction-Opening Curves Method

Traction-Opening Curves can be directly defined on the MODEL field of the MCOHE Bulk Data Entry.

CIFHEX and CIFPEN entries are used for modeling cohesive elements. The PCOHE property defines properties of cohesive elements. The thickness of cohesive element layer can be non-zero in the actual geometry, but is internally always considered to be 1.0.

For the Traction-Opening Curves Method, only a single layer of cohesive elements are allowed to model a particular adhesive or bonded interface.

Typically, tensile and shear deformation at the interface are of interest to determine the integrity and/or degradation of the adhesion/bonding. Stiffness in compression is controlled via the SFC field on the MCOHE entry.

The relative displacement between the nodes of the top and bottom faces are calculated.

The three displacements (

d_{x}

d_{y}

, and

d_{z}

) are then combined. The normal and shear modes are mixed based on displacement.(1)

d_{e f f} = \sqrt{{(β d_{x})}^{2} + {(β d_{y})}^{2} + {(\max \{0.0, d_{z}\})}^{2}}

Where,

$d_{x}$ , $d_{y}$ , and $d_{z}$: Are the relative displacements of the top and the bottom faces of a cohesive element along elemental x-, y-, and z-axes.
$β$: Is the mixing coefficient, which can be input on the BETA field.

Using this combined relative displacement ( $d_{e f f}$ ), the combined traction is determined based on the chosen traction-opening curve (MODEL filed on MOCHE).

Damage Model Method

Damage Model method allows modeling of finite thickness of adhesive/bonded interfaces. This allows for modeling adhesive surfaces based on experimental data. In addition, damage initiation and propagation/evolution can be modeled.

MCOHED, DMGINI (Damage Initiation), and DMGEVO (Damage Evolution) entries are mandatory. Depending on the modeling technique (Element-based vs Contact-based), either CIFHEX/CIFPEN/PCOHE entries or CONTACT interfaces are required. Refer to Modeling Techniques.

Penalty stiffnesses in three directions can be defined on the KZ, KX, and KY fields of the MCOHED entry, in which KZ is for normal direction, KX and KY are for two tangential directions. Single or multiple layers of cohesive elements can be defined in the interface.

For Element-based modeling, the thickness of the cohesive element layer can be defined using the THICKNESS field on PCOHE. For Contact-based modeling, the thickness of the cohesive zone is considered internally to be equal to 1.0.

The DMGINID and DMGEVOID fields on MCOHED can be used to identify the mandatory DMGINI and DMGEVO Bulk Data Entries in element-based modeling. When contact based modeling is used, stiffness in compression is determined by contact property.

The relative displacement between the nodes of the top and bottom faces are calculated (similar to the traction-opening method). First, the trial traction values ( $K_{i} d_{i}$ ) are calculated by multiplying the elastic moduli and relative displacements in each of the three directions.

Next, the determination of damage initiation is carried out, using the specified criteria on the CRI field of the DMGINI entry.

Strain-Based Criteria

The maximum strain value is defined on the V1, V2, V3 fields of the DMGINI entry.
The actual strain is calculated by the formula: (relative displacement divided by thickness). Where, the thickness is defined by the THICKNESS field on PCOHE.
Using both maximum strain and actual strain, the Damage initiation determination is done based on the following formula:

MAXE

(2) $m a x {\frac{ε_{x}}{m a x e_{x}}, \frac{ε_{y}}{m a x e_{Y}}, \frac{ε_{z}}{m a x e_{Z}}} = 1$
Where,

$m a x e_{x}$

Is V1.

$m a x e_{Y}$

Is V2.

$m a x e_{Z}$

Is V3.

QUADE

(3) ${(\frac{ε_{x}}{m a x e_{x}})}^{2} + {(\frac{ε_{y}}{m a x e_{Y}})}^{2} + {(\frac{ε_{z}}{m a x e_{Z}})}^{2} = 1$
Where,

$m a x e_{x}$

Is V1.

$m a x e_{Y}$

Is V2.

$m a x e_{Z}$

Is V3.

(4) $ε_{i} = \frac{d_{i}}{t_{0}}$
Where, $i = x, y, z$ and $t_{0}$ are the thickness defined in the PCOHE entry.

Stress-Based Criteria

The maximum stress value is defined on the V1, V2, V3 fields of the DMGINI entry.
The actual stress is the value of trial traction in each of the corresponding three directions.
Using both maximum stress and actual stress, the damage initiation determination is done based on the following formula:

MAXS

(5) $m a x {\frac{σ_{x}}{m a x σ_{x}}, \frac{σ_{y}}{m a x σ_{Y}}, \frac{σ_{z}}{m a x σ_{Z}}} = 1$
Where,

$m a x σ_{x}$

Is V1.

$m a x σ_{Y}$

Is V2.

$m a x σ_{Z}$

Is V3.

(6) $σ_{i} = k_{i} \frac{d_{i}}{t_{0}}$
Where,

$i = x, y, z$ $t_{0}$

Are the thickness defined in the PCOHE entry.

$k_{i}$

Is the initial elastic stiffness.

QUADS

(7) ${(\frac{σ_{x}}{m a x σ_{x}})}^{2} + {(\frac{σ_{y}}{m a x σ_{Y}})}^{2} + {(\frac{σ_{z}}{m a x σ_{Z}})}^{2} = 1$
Where,

$m a x σ_{x}$

Is V1.

$m a x σ_{Y}$

Is V2.

$m a x σ_{Z}$

Is V3.

In the formula above, the x-, y- and z-directions are in material coordinate system (refer to MCOHED card). The z-direction is the normal direction and x- and y- directions are two tangential directions.

If damage initiation criteria are not satisfied, there is no damage. Trial traction is equal to actual traction.

Therefore, there is no crack initiation and propagation and the corresponding cohesive-related output are printed to the result files.

If damage initiation criteria are satisfied, this implies that damage is initiated. The damage index is now calculated. There are two methods of damage index calculation:

Displacement-Based Damage Index (TYPE=COHDISP on DMGEVO entry)
Energy Dissipation-Based Damage Index (TYPE=COHENRG on DMGEVO entry)

For either type of damage index calculation, either linear (SHAPE=LIN) or exponential (SHAPE=EXP), shapes of the traction-opening curve can be used on the DMGEVO entry.

Displacement-Based Damage Index

If SHAPE = LIN:

If the traction

T_{i}

decreases linearly during damage evolution, then the damage index is formulated by:(8)

D = \frac{d_{f}}{d_{m a x}} \frac{(d_{m a x} - d_{o})}{(d_{f} - d_{o})}

If SHAPE=EXP

If the traction

T_{i}

decreases exponentially during damage evolution, then the damage index is formulated by:(9)

D = 1 - \frac{d_{o}}{d_{m a x}} [1 - \frac{1 - exp (- α \frac{d_{m a x} - d_{o}}{d_{f} - d_{o}})}{1 - exp (- α)}]

Where,

$d_{m a x}$: Is the maximum separation ( $d$ ) in history. If the cohesive zone is only being loaded, then $d_{m a x}$ is equal to the current $d$ , which is calculated by OptiStruct and updated at each step. If cohesive zone is also unloaded, then within the unloading zone, $d_{m a x}$ is equal to the maximum $d$ in history (as it could be a previous value of $d$ , prior to the beginning of unloading).
$d_{o}$: Is the critical separation (the separation, $d$ , when damage is initiated – when the crack initiation criteria are satisfied).
$d_{f}$: Is the maximum separation ( $d_{f} = d_{o} + W_{1}$ ).
$α$: Is the ALPHA field on DMGEVO entry.
$d$: Is the current opening at each step of the solution.; (10) $d = \sqrt{{(d_{x})}^{2} + {(d_{y})}^{2} + m a x {(0.0, d_{z})}^{2}}$
$W_{1}$: Is the W1 field on the DMGEVO entry.
$d_{x}$ , $d_{y}$ , and $d_{z}$: Are the relative displacements of the top and the bottom faces of a cohesive element along elemental x-, y-, and z-axes.

Energy Dissipation-Based Damage Index

For Energy Dissipation-based Damage Index, the critical total energy ( $G_{c}$ ) is the key value used for calculations. Its calculation and usage depend on the type of curve (LIN/EXP) and the mode mix method (MMXFM = blank, 1,2). $G_{c}$ represents the energy at which failure occurs.

W1, W2, and W3 fields on the DMGEVO entry define the critical energies in each of the three fracture modes (normal, in-plane shear, and out-of-plane shear, respectively).

The following nomenclature help identify the fracture modes vs the corresponding directions.

Normal Fracture mode:: Mode 1; Direction Z; W1 = $G_{I c}$ (power law) = $G_{n c}$ (Benzeggah-Kenane (BK) Form) defines the critical energy for normal failure. This is the energy at which failure occurs in normal direction.; $G_{I}$ (power law) = $G_{n}$ (BK form) defines the energy (area under the curve) until the current step of the solution in the normal direction.
In-Plane Shear Fracture mode:: Mode 2; Direction X; W2 = $G_{I I c}$ (power law) defines the critical energy for in-plane shear failure. This is the energy at which failure occurs in in-plane-shear direction.; W2 = $G_{s c}$ (BK Form) defines the critical energy for shear failure (it is interpreted as a combined resultant shear for BK Form). This is the energy at which failure occurs in the resultant shear direction.; $G_{I I}$ (Power law) defines the energy (area under the curve) until the current step of the solution in the in-plane shear direction.; $G_{s}$ (BK form) defines the energy (area under the curve) until the current step of the solution in the shear direction (it is interpreted as a combined resultant shear for BK form).
Out-of-Plane Shear Fracture mode:: Mode 3; Direction Y; W3 = $G_{I I I c}$ (power law) defines the critical energy for out-of-plane shear failure. This is the energy at which failure occurs in out-of-plane shear direction.; $G_{I I I}$ defines the energy (area under the curve) until the current step of the solution in the out-of-plane shear direction.

Note: For BK Form, both in-plane and out-of-plane shear are combined and are considered as a resultant shear. W3 is not applicable for BK form. So

G_{s}

is considered as the energy in resultant shear for BK form (unlike separate

G_{I I}

energy in in-plane shear and

G_{I I I}

energy in out-of-plane shear for Power law), while

G_{n}

is the energy in normal direction for BK form, similar to

G_{I}

, which is the energy in normal direction for Power Law.

$G_{o}$ is the elastic energy absorbed by cohesion. It is the area under the straight line section (before damage initiation) of the exponential curve.

The compound traction is:(11)

T = \sqrt{T_{x}^{2} + T_{y}^{2} + T_{z}^{2}}

Where, $D$ is the damage index (it is always ≤ 1.0).

If SHAPE = LIN:

If the traction

T_{i}

decreases linearly during damage evolution, then the damage index is formulated by:(12)

D = \frac{d_{f}}{d_{m a x}} \frac{(d_{m a x} - d_{o})}{(d_{f} - d_{o})}

Where,

$d_{m a x}$: Is the maximum separation ( $d$ ) in history. If the cohesive zone is only being loaded, then $d_{m a x}$ is equal to the current $d$ , which is calculated by OptiStruct and updated at each step. If cohesive zone is also unloaded, then within the unloading zone, $d_{m a x}$ is equal to the maximum $d$ in history (as it could be a previous value of $d$ , prior to the beginning of unloading).
$d_{o}$: Is the critical separation (the separation, $d$ , when damage is initiated – when the crack initiation criteria are satisfied).
$d_{f}$: Is the separation at which zero traction is produced in analysis.; Is estimated by the released energy ( $G_{c}$ ):(13) $d_{f} = \frac{2 G_{c}}{T_{\max}}$
$T_{\max}$: Is the effective traction when crack initiation criteria are satisfied.
$G_{c}$: Depends on the mode mix form (MMXFM field on DMGEVO entry).

If MMXFM field is blank:(14) $G_{c} = W_{1}$
If MMXFM field is set to 1 (Power Law):
The $G_{c}$ is derived from:(15) ${(\frac{G_{I}}{G_{I c}})}^{α} + {(\frac{G_{I I}}{G_{I I c}})}^{α} + {(\frac{G_{I I I}}{G_{I I I c}})}^{α} = {(\frac{G}{G_{c}})}^{α}$ Where, $G_{I c}$ , $G_{I I c}$ , and $G_{I I I c}$ are the W1, W2, and W3 fields on the DMGEVO entry.

The value of $G_{c}$ is derived as:(16) $G_{c} = d^{2} {({(\frac{d_{x}^{2}}{W_{2}})}^{α} + {(\frac{d_{y}^{2}}{W_{3}})}^{α} + {(\frac{\max {(0.0, d_{z})}^{2}}{W_{1}})}^{α})}^{- \frac{1}{α}}$

Where,

$α$

Is the ALPHA field on DMGEVO entry.

$W_{1}$ , $W_{2}$ , and $W_{3}$

Are the W1, W2, and W3 fields on the DMGEVO entry.

$d_{x}$ , $d_{y}$ , and $d_{z}$

Are the relative displacements of the top and the bottom faces of a cohesive element along elemental x-, y-, and z-axes.
If MMXFM field is set to 2 (BK Form):
The value of $G_{c}$ is based on the following equation:(17) $G_{n c} + (G_{s c} - G_{n c}) {(\frac{G_{s}}{G_{T}})}^{η} = G_{c}$

Where,

$η$

Is the ALPHA field on DMGEVO entry.

$G_{n c}$ and $G_{s c}$

Are the W1 and W2 fields on the DMGEVO entry.

If SHAPE=EXP:

If the traction

T_{i}

decreases exponentially during damage evolution, then the damage index is formulated by:(18)

D = \int_{d_{o}}^{d_{f}} \frac{T}{G_{c} - G_{o}} d d

Where,

$T$: Is the compound traction.
$G_{o}$: Is the elastic energy absorbed by the cohesion.
$d_{o}$: Is the critical separation (the separation, $d$ , when damage is initiated – when the crack initiation criteria are satisfied).
$d_{f}$: Is the final separation.
$G_{c}$: Is the total energy that can be dissipated by cohesion under the current opening pattern (combination of $d_{x}$ , $d_{y}$ , and $d_{z}$ ). The released energy, $G_{c}$ , is calculated automatically by OptiStruct and depends on the mode mix form (MMXFM field on DMGEVO entry) and the energy that can be dissipated in each mode (W1, W2, and W3).

If MMXFM field is blank:(19) $G_{c} = W_{1}$
If MMXFM field is set to 1 (Power Law):
The $G_{c}$ is derived from:(20) ${(\frac{G_{I}}{G_{I c}})}^{α} + {(\frac{G_{I I}}{G_{I I c}})}^{α} + {(\frac{G_{I I I}}{G_{I I I c}})}^{α} = {(\frac{G}{G_{c}})}^{α}$

Where,

$G_{I c}$ , $G_{I I c}$ , and $G_{I I I c}$

Are the W1, W2, and W3 fields on the DMGEVO entry.

$G_{I}$ , $G_{I I}$ , and $G_{I I I}$

Are the energies under the curve until the current step. They depend on the type of curve (LIN/EXP).

Therefore, the derivation of $G_{c}$ also depends on the type of curve.

However, by default, the same value of $G_{c}$ as the linear curve is used for the exponential curve for Power law.(21) $G_{c} = d^{2} {({(\frac{d_{x}^{2}}{W_{2}})}^{α} + {(\frac{d_{y}^{2}}{W_{3}})}^{α} + {(\frac{\max {(0.0, d_{z})}^{2}}{W_{1}})}^{α})}^{- \frac{1}{α}}$

Where,

$α$

Is the ALPHA field on DMGEVO entry.

$W_{1}$ , $W_{2}$ , and $W_{3}$

Are the W1, W2, and W3 fields on the DMGEVO entry.

$d_{x}$ , $d_{y}$ , and $d_{z}$

Are the relative displacements of the top and the bottom faces of a cohesive element/contact along elemental/local x-, y-, and z-axes.
If MMXFM field is set to 2 (BK Form):
The value of $G_{c}$ is:(22) $G_{n c} + (G_{s c} - G_{n c}) {(\frac{G_{s}}{G_{T}})}^{η} = G_{c}$

Where,

$η$

Is the ALPHA field on DMGEVO entry.

$G_{n c}$ and $G_{s c}$

Are the W1 and W2 fields on the DMGEVO entry.

Calculating Actual Traction

Actual traction is calculated as follows, based on the Damage Index calculation selected above.

When DMGEVO is referenced by MCOHED, the traction in cohesive elements is calculated based on the damage index,

D

.(23)

T_{i} = (1 - D) k_{i} d_{i}

(24)

i = x, y, z

Where,

$k_{i}$ and $d_{i}$: Are the initial elastic stiffness and opening in the i-direction respectively.
$t_{0}$: Is the thickness defined in the PCOHE entry.

This actual traction is then subsequently used for the solution.

Cohesive Element Erosion

When damage index of all integration points in a cohesive element reaches the value defined on the MXDMG field in MCOHED card, and none of these integration points is in compression, the cohesive element is eroded and do not work in the rest analysis in the current subcase and the continued subcases.

Modeling Techniques

There are currently two techniques to activate adhesion/bonding and simulate the cohesive behavior at the corresponding interfaces:

Defining Cohesive Elements Directly (CIFHEX/CIFPEN elements)
Contact-based method (CONTACT entry)

Defining Cohesive Elements (CIFHEX and CIFPEN elements)

The definition of cohesive elements (CIFHEX and CIFPEN) can be used For the Traction-Opening Curves method and the Damage Models method. The formulation of these two cohesive elements are as follows:

The main focus of the CIFHEX and CIFPEN elements is the relative movement between the top and bottom faces.
The relative displacement between the nodes of the top and bottom faces at each integration point in each of the three directions (elemental X, Y, Z) determines the cohesive opening.

Figure 5.
Traction provides the tensile and shear stiffness and SFC field on the PCOHE entry identifies the compressive stiffness for cohesive elements.
Note: Traction calculation depends on the method used for Adhesive/Bonded interface modeling.
Refer to Interface Elements for more information on the cohesive element formulation.

The recommendations to define cohesive elements:

Cohesive elements should be inserted in the path of crack propagation.
For Traction-Opening Method, only a single layer of cohesive elements should be used.
For Damage Model Method, multiple layers of cohesive elements can be used.
CIFHEX and CIFPEN elements are available for modeling cohesive elements.
Cohesive elements can only be connected to shell or solid elements of the base model.
If an exact one-on-one nodal correspondence, along with the exact same mesh density exists between the cohesive element layer and the shell/solid connecting layer of the base model, then the nodes can be shared (equivalenced) and no contact definition is required.
If such an exact one-on-one nodal correspondence does not exist, the CONTACT(FREEZE) or TIE connection should be used to connect the cohesive elements on either top/bottom layer to the corresponding shell/solid elements of the base model.
Cohesive elements may have a geometrical thickness in the interface. For the Traction-Opening method, a thickness of 1.0 is internally used automatically, regardless of the geometric thickness. For the Damage Model method, the THICKNESS field on PCOHE can be used to control the thickness interpretation.
In some cases, it may be difficult for analysis with cohesive elements to converge. Damping stabilization can be introduced in the cohesive elements to help convergence. The damping stabilization can be defined by VED on MCOHE and MCOHED entries. Damping stabilization is currently not available in cohesive contact.

Contact-Based Technique (CONTACT entry)

The Contact-based technique does not require the use of cohesive elements (CIFHEX/CIFPEN) to model the cohesive zone. This technique allows simplification of the model setup, as it eliminates the need to mesh and setup the cohesive elements.

The COHE continuation line on the CONTACT Bulk Data Entry can be used to activate the Contact-based method for cohesive zone modeling. The MCOHEDID field references the MCOHED identification number and therefore identifies the contact interface as an adhesive/bonded interface.

Only Damage Model Method is available for CONTACT-based modeling. Also, the thickness of the cohesive zone is internally always set to 1.0 for this method.

Contact penalty is used to avoid penetration in compression. Contact effects are ignored, and cohesive effects are activated in separation.

The COHEGSET field on the COHE continuation line defines the initial configuration of the contact interface.

Note: This is only applied to the initial state of the analysis and as the analysis progresses, the contact definition varies normally as with any nonlinear analysis.

Currently only SMALL sliding, frictionless, N2S/S2S contact is supported for cohesive modeling.

Supported Solution Sequences

Cohesive elements are currently supported for the following solutions:

Nonlinear Static Analysis (SMDISP and LGDISP)
Nonlinear Transient Analysis (SMDISP and LGDISP)
- Mass is not considered in cohesive elements.
Linear Analysis, including static, transient, buckling, and eigen mode analysis
- Cohesive effects are not available for Linear Analysis.
- The initial stiffness of cohesive elements is used in linear analysis. The initial stiffness is determined by the MCOHE or MCOHED entry (the initial slope of the traction-separation curve defined in MCOHE entry or the Ki value defined in the MCOHED entry.
- There is no crack propagation/initiation.
- There are no cohesive elements related output for Linear Analysis.
Cohesive elements are currently only supported for Implicit Analysis. Explicit Analysis is not supported.

Output

Output from the cohesive zone is currently only available in H3D format.

The following results are specifically output for cohesive zones:

Cohesive Damage Initiation Index
- This item is shown in cohesive elements and secondary surface of cohesive contact.
Cohesive Damage Index (Damage Index)
- This item is shown in cohesive elements and secondary surface of cohesive contact.
Cohesive Energy by Mode (Dissipated Cohesive Energy by Mode)
- Cohesive Energy is output in terms of 3 modes (Mode I, II, and III)
- This item is not available in cohesive contact
Cohesive Energy per Area by Mode (Dissipated Cohesive Energy per Area by Mode)
- Cohesive Energy per area is output in terms of 3 modes (Mode I, II, and III)
- This item is not available in cohesive contact
Cohesive Maximum Opening in History (Maximum Opening)
- This is the maximum relative displacement in history.
- This item is not available in cohesive contact.
Cohesive Opening(s)
- This provides output for relative displacement in local elemental system and basic system.
Cohesive Status (Status)
- Indicates element loading/unloading/fail state
- 0: Loading
- 1: Unloading/Reloading
- 2: Fail
Cohesive Traction (Traction)

Eroded cohesive elements are not shown in h3d output from its eroded analysis time. Therefore, the above items are not available in eroded cohesive elements.

In cohesive zone output when contacts are used to model the cohesive zone, cohesive traction and cohesive opening are listed in label ‘Contact Traction / Normal’, ‘Contact Traction / Tangent’, ‘Contact Deformation / Normal’ and ‘Contact Deformation/Tangent.’ In order to keep it consistent with contact pressure, cohesive normal traction is shown as negative value in ‘Contact Traction / Normal.’