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/techniques by which interfaces can currently be modeled in OptiStruct. In the damagebased method that relies on cohesive elements, thickness of the adhesive/bonded interface can be considered. For other methods/techniques, thickness effect is excluded. The term technique is used to distinguish between the element types that are used, which are cohesive elements and contact elements.
Interface Implementation Methods
Cohesive zones are used to model adhesive/bonded interfaces, where layers/parts are connected with resin or adhesive.
 Potentialbased method
 Damagebased method
There are three modes for the cohesive zone deformation.
The system xyz in Figure 2 is the cohesive zone coordinate system. Unless it is specified in the PCOHE card, the establishment is shown in CIFHEX/CIFPEN Bulk Data Enty. Mode I is the normal opening in zdirection, Mode II and III are the shear opening in xz plane and yz plane, respectively.
In contact, it is assumed that contact pairs have the same property in two tangential directions. Therefore, it is suggeted that the cohesive property in Mode II and III be the same.
Potentialbased method (MCOHE)
In the potentialbased method the tractionopening model is selected between three types of curves.
In the potentialbased method, only a single layer of cohesive elements is allowed for modeling a particular adhesive or bonded interface.
Typically, tensile and shear deformations 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}_{I},{d}_{II}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{d}_{III}$ ) are used to derive the combined relative displacement with the mixing formulation:
 ${d}_{I},{d}_{II}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{d}_{III}$
 The opening in Mode I, II and III.
 $\beta $
 The mixing coefficient, which can be input on the BETA field.
 Bilinear:$$T=\{\begin{array}{cc}\frac{2G}{{d}_{m}\text{}}\frac{d}{{d}_{c}}& 0\le d\le {d}_{c}\\ \frac{2G}{{d}_{m}}\left(\frac{{d}_{m}d}{{d}_{m}{d}_{c}}\right)& {d}_{c}\le d\le {d}_{m}\\ 0& d{d}_{m}\end{array}$$
 Exponential:$$T=G\frac{d}{{d}_{c}^{2}}{e}^{\frac{d}{{d}_{c}}}$$
 Linearexponential:$$T=\{\begin{array}{cc}\frac{2qG}{{d}_{c}\text{(q+2)}}\frac{d}{{d}_{c}}& 0\le d\le {d}_{c}\\ \frac{2qG}{{d}_{c}\text{(q+2)}}{e}^{q\left(1\frac{d}{{d}_{c}}\right)}& d{d}_{c}\end{array}$$
 ${d}_{c}$
 is the CRTOD
 ${d}_{m}$
 is the MAXOD
 $d$
 is the combined relative displacement ( ${d}_{eff}$ )
Also, in the linearexponential curve definition, one more parameter, $q$ , is required, which is input in the field EXP in MCOHE card.
Damagebased Method (MCOHED)
Damagebased method allows modeling of finite thickness of adhesive/bonded interfaces. This allows for modeling adhesive layers based on experimental data.
In this method, crack growth is controlled by development of damage in each element. The damage of an element is controlled by two indices, namely, damage initiation index and damage evolution index. In each element, the two indices are initially zero. With loading, damage initiation index is increasing. When damage initiation index reaches 1.0, damage appears. Damage initiation index keeps being 1.0 afterwards and damage evolution index starts to increase. When damage evolution index reaches 1.0, the damage is matured and there is no cohesion afterwards. Thus, crack advances.
The damage initiation and evolution are controlled by several parameters defined in MCOHED, DMGINI (Damage Initiation) and DMGEVO (Damage Evolution) cards by users.
Depending on the modeling technique (elementbased versus contactbased), either CIFHEX/CIFPEN/PCOHE entries or CONTACT interfaces are required. Refer to Modeling Techniques.
Elasticity modulus/Penalty stiffness in the three directions are defined on the KI, KII, and KIII fields of the MCOHED entry, in which KI is for the normal direction, KII and KIII are for the two tangential directions. Single or multiple layers of cohesive elements can be defined in the interface.
Typically, tensile and shear deformations 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 MCOHED entry in the elementbased technique. When contactbased technique is used, stiffness in compression is determined by contact property.
For the elementbased technique, thickness of the cohesive element layer can be defined using the THICKNESS field on PCOHE. For contactbased technique, thickness of the cohesive zone is considered internally to be one unit.
The following introduction on damagebased cohesive model is for elementbased technique. They are also valid for contactbased technique, if thickness $t$ is replaced by unity.
The DMGINID and DMGEVOID fields on MCOHED are used to specify the mandatory DMGINI and DMGEVO Bulk Data Entries in elementbased 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 potentialbased method. First, the trial traction values ( ${K}_{i}{d}_{i}$ / ${t}_{0}$ ) are calculated by multiplying the penalty stiffness (which is elasticity modulus divided by thickness in quantity) and the opening in the three modes.
Next, determination of the damage initiation is carried out, using the specified criteria on the CRI field of the DMGINI entry.
Strainbased Initiation 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 the thickness)
Where, 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

$$max\left\{\frac{{\epsilon}_{I}}{max\left({e}_{I}\right)},\frac{{\epsilon}_{II}}{max\left({e}_{II}\right)},\frac{{\epsilon}_{III}}{max\left({e}_{III}\right)}\right\}=1$$Where,
 $max{e}_{\text{I}}$
 = V1.
 $max{e}_{\mathrm{II}}$
 = V2.
 $max{e}_{\mathrm{III}}$
 = V3.
 QUADE

$${\left(\frac{{\epsilon}_{I}}{max\left({e}_{I}\right)}\right)}^{2}+{\left(\frac{{\epsilon}_{II}}{max\left({e}_{II}\right)}\right)}^{2}+{\left(\frac{{\epsilon}_{III}}{max\left({e}_{III}\right)}\right)}^{2}=1$$Where,
 $max{e}_{\text{I}}$
 = V1.
 $max{e}_{\mathrm{II}}$
 = V2.
 $max{e}_{\mathrm{III}}$
 = V3.
Stressbased Initiation 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

$$max\left\{\frac{{\sigma}_{I}}{max\left({\sigma}_{I}\right)},\frac{{\sigma}_{II}}{max\left({\sigma}_{II}\right)},\frac{{\sigma}_{III}}{max\left({\sigma}_{III}\right)}\right\}=1$$Where,
 $max{\sigma}_{\text{I}}$
 = V1.
 $max{\sigma}_{\mathrm{II}}$
 = V2.
 $max{\sigma}_{\mathrm{III}}$
 = V3.
 QUADS

$${\left(\frac{{\sigma}_{I}}{max\left({\sigma}_{I}\right)}\right)}^{2}\text{}+{\left(\frac{{\sigma}_{II}}{max\left({\sigma}_{II}\right)}\right)}^{2}+{\left(\frac{{\sigma}_{III}}{max\left({\sigma}_{III}\right)}\right)}^{2}=1$$Where,
 $max{\sigma}_{I}$
 = V1.
 $max{\sigma}_{\mathrm{II}}$
 = V2.
 $max{\sigma}_{\mathrm{III}}$
 = V3.
If damage initiation criteria are not satisfied, then there is no damage. The trial traction is equal to actual traction. Therefore, there is no crack initiation and propagation, and the corresponding cohesiverelated output are printed to the result files.
 Displacementbased damage evolution index (TYPE=COHDISP on DMGEVO entry)
 Energy dissipationbased damage evolution index (TYPE=COHENRG on DMGEVO entry)
For either type of damage evoluation index calculation, either linear (SHAPE=LIN) or exponential (SHAPE=EXP), shapes of the tractionopening curve can be used on the DMGEVO entry.
$D$ is the damage evolution index (it is always ≤ 1.0) and in the output it is referred as damage index.
Displacementbased Damage Evolution Index
If SHAPE = LIN:
If the traction ${T}_{i}$ decreases linearly during damage evolution, then the damage index is formulated by:
If SHAPE=EXP
 ${d}_{max}$
 Maximum opening ( $d$ ) in history. If the cohesive zone is only being loaded, then ${d}_{max}$ 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}_{max}$ 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}$
 Critical opening (the opening, $d$ , when damage is initiated – when the crack initiation criteria are satisfied).
 ${d}_{f}$
 Maximum opening ( ${d}_{f}={d}_{o}+{W}_{1}$ ).
 ${W}_{1}$
 W1 field on the DMGEVO entry.
 $\alpha $
 ALPHA field on DMGEVO entry.
 $d$
 Current opening at each step of the solution.
Energy Dissipationbased Damage Evolution Index
For energy dissipationbased damage evoluation index, the critical total energy ( ${G}_{c}$ ) also known as fracture toughness is the key value used for calculations. Its calculation and usage depend on the type of the 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 Fracture Mode I:
 W1 = ${G}_{Ic}$ (Power law) = ${G}_{nc}$ (BenzeggahKenane (BK) form) defines the critical energy at which failure occurs in normal direction.
 Shear Fracture Mode II:
 W2 = ${G}_{IIc}$ (Power law) defines the critical energy at which failure occurs in Mode II.
 Shear Fracture Mode III:
 W3 = ${G}_{IIIc}$ (Power law) defines the critical energy at which failure occurs in Mode III.
${G}_{T}={G}_{n}+{G}_{s}$ is the total energy (area under the curve) until the current solution step.
If SHAPE = LIN:
If the traction ${T}_{i}$ decreases linearly during the damage evolution, then the damage evolution index is formulated by:
 ${d}_{max}$
 The maximum opening ( $d$ ) in history. If the cohesive zone is only being loaded, then ${d}_{max}$ 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}_{max}$ 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}$
 The critical opening (the opening, $d$ , when damage is initiated – when the crack initiation criteria are satisfied).
 ${d}_{f}$
 The opening at which zero traction is produced in the analysis.
 ${T}_{\mathrm{max}}$
 Effective traction when crack initiation criteria are satisfied.
 If MMXFM field is blank:$${G}_{c}={W}_{1}$$
 If MMXFM field is set to 1 (Power law):
The ${G}_{c}$ is derived from:
$${\left(\frac{{G}_{I}}{{G}_{Ic}}\right)}^{\alpha}+{\left(\frac{{G}_{II}}{{G}_{IIc}}\right)}^{\alpha}+{\left(\frac{{G}_{III}}{{G}_{IIIc}}\right)}^{\alpha}={\left(\frac{G}{{G}_{c}}\right)}^{\alpha}$$Where, ${G}_{Ic}$ , ${G}_{IIc}$ , and ${G}_{IIIc}$ are the W1, W2, and W3 fields on the DMGEVO entry.The value of ${G}_{c}$ is derived as:
$${G}_{c}={\left[{\left(\frac{{d}_{II}}{d}\right)}^{2\alpha}{\left(\frac{1}{{W}_{2}}\right)}^{\alpha}+{\left(\frac{{d}_{III}}{d}\right)}^{2\alpha}{\left(\frac{1}{{W}_{3}}\right)}^{\alpha}+{\left(\frac{\mathrm{max}\left(0.0,{d}_{I}\right)}{d}\right)}^{2\alpha}{\left(\frac{1}{{W}_{I}}\right)}^{\alpha}\right]}^{\frac{1}{\alpha}}$$Where, $\alpha $ is the ALPHA field on DMGEVO entry.
 If MMXFM field is set to 2 (BK form):
The value of ${G}_{c}$ is based on the following equation:
$${G}_{nc}+\left({G}_{sc}{G}_{nc}\right){\left(\frac{{G}_{s}}{{G}_{T}}\right)}^{\eta}={G}_{c}$$Where, $\eta $
 ALPHA field on DMGEVO entry.
 ${G}_{nc}$ and ${G}_{sc}$
 W1 and W2 fields on the DMGEVO entry.
${G}_{T}={G}_{n}+{G}_{s}$ is the total energy (area under the curve) until the current solution step.
If SHAPE=EXP:
If the traction ${T}_{i}$ decreases exponentially during the damage evolution, then the damage evolution index is formulated by:
 $T$
 Compound traction.
 ${G}_{o}$
 Elastic energy absorbed by the cohesion. It is the area under the straight linesection (before damage initiation) of the exponential curve.
 ${d}_{o}$
 Critical opening (the opening, $d$ , when damage is initiated – when the crack initiation criteria are satisfied).
 ${d}_{f}$
 Final opening.
 ${G}_{c}$
 Total energy that can be dissipated by cohesion under the current opening pattern (combination of ${d}_{\text{I}}$ , ${d}_{\mathrm{II}}$ , and ${d}_{\mathrm{III}}$ ). The critical energy, ${G}_{c}$ , is calculated automatically by OptiStruct and it 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:$${G}_{c}={W}_{1}$$
 If MMXFM field is set to 1 (Power law):
The ${G}_{c}$ is derived from:
$${\left(\frac{{G}_{I}}{{G}_{Ic}}\right)}^{\alpha}+{\left(\frac{{G}_{II}}{{G}_{IIc}}\right)}^{\alpha}+{\left(\frac{{G}_{III}}{{G}_{IIIc}}\right)}^{\alpha}={\left(\frac{G}{{G}_{c}}\right)}^{\alpha}$$Where, ${G}_{Ic}$ , ${G}_{IIc}$ , and ${G}_{IIIc}$
 The W1, W2, and W3 fields on the DMGEVO entry.
 ${G}_{I}$ , ${G}_{II}$ , and ${G}_{III}$
 The energies under the tractionopening curve until the current step. They depend on the type of curve (LIN/EXP).
However, by default, for the exponential curve, the same value of ${G}_{c}$ is used for the linear curve.
$${G}_{c}={\left[{\left(\frac{{d}_{II}}{d}\right)}^{2\alpha}{\left(\frac{1}{{W}_{2}}\right)}^{\alpha}+{\left(\frac{{d}_{III}}{d}\right)}^{\alpha}{\left(\frac{1}{{W}_{3}}\right)}^{\alpha}+{\left(\frac{\mathrm{max}\left(0.0,{d}_{I}\right)}{d}\right)}^{2\alpha}{\left(\frac{1}{{W}_{I}}\right)}^{\alpha}\right]}^{\frac{1}{\alpha}}$$  If MMXFM field is set to 2 (BK form):
The value of ${G}_{c}$ is:
$${G}_{nc}+\left({G}_{sc}{G}_{nc}\right){\left(\frac{{G}_{s}}{{G}_{T}}\right)}^{\eta}={G}_{c}$$Where, $\eta $
 ALPHA field on DMGEVO entry.
 ${G}_{nc}$ and ${G}_{sc}$
 W1 and W2 fields on the DMGEVO entry.
${G}_{T}={G}_{n}+{G}_{s}$ is the total energy (area under the curve) until the current solution step.
Calculating Actual Traction
Actual traction is calculated as follows, based on the damage evolution index calculation described in the previous section.
When DMGEVO is referenced by MCOHED, the traction in cohesive elements is calculated based on the damage evolution index, $D$ .
$i=\text{I},\mathrm{II},\mathrm{III}$
 ${k}_{i}$ and ${d}_{i}$
 The elasticity modulus and opening corresponding to Mode $i$ .
 ${t}_{0}$
 Thickness defined in the PCOHE entry.
This actual traction is then subsequently used for the solution.
Cohesive Element Erosion
When damage evoluation 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 does not work in the rest of the analysis in the current subcase and the continued subcases.
Additionally, eroded elements are not shown in H3D output starting from when they are eroded during the analysis.
Modeling Techniques
 Elementbased technique (CIFHEX/CIFPEN elements)
 Contactbased technique (CONTACT entry)
Elementbased Technique (CIFHEX and CIFPEN elements)
 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.
 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.
 Cohesive elements should be inserted in the path of crack propagation.
 For potentialbased method, only a single layer of cohesive elements should be used.
 For damagebased 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 onetoone nodal correspondence exists along with exactly the same mesh density 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 onetoone 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 potentialbased method, unit thickness is internally used automatically, regardless of the geometric thickness and thus, the thickness effect is excluded. For the damagebased method, the THICKNESS field on PCOHE can be used to control the thickness interpretation. It is important to use realistic value. The default value is set to 1.0.
 In some cases, it may be difficult to obtain convergence with cohesive elements. 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.
ContactBased Technique (CONTACT entry)
The contactbased 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 contactbased 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 damagebased method is available for contactbased technique. Also, the thickness of the cohesive zone is internally always set to unit for this technique. Thus, the thickness effect is excluded.
Contact penalty is used to avoid penetration in compression in opening cohesive effects are activated and contact effects are ignored
Currently only SMALL sliding, frictionless, N2S/S2S contact is supported for cohesive modeling.
Supported Solution Sequences
 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 tractionopening 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.
 Cohesive Damage Initiation Index
 Damage initiation index is calculated at the center of a cohesive element. The damage initiation index for each cohesive element indicates how far away the current state of an element is from its peak traction.
 Damage for an element is considered to have begun when the Cohesive damage initiation index value reaches 1.0. At this point the Cohesive damage index starts to grow above 0.0.
 This item is shown in cohesive elements and secondary surface of the cohesive contact.
 Cohesive Damage Index (Damage Evolution Index)
 This value is calculated for each cohesive element and it starts to grow above 0.0 when the damage for the element begins (that is, when the cohesive damage initiation index reaches 1.0).
 This item is shown in cohesive elements and secondary surface of the cohesive contact.
 Cohesive Energy per Area by Mode (Dissipated Cohesive Energy per Area by
Mode)
 This is the area under the traction opening curve for each 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 Energy by Mode (Dissipated Cohesive Energy by Mode)
 This is the “Cohesive energy per area by mode” multiplied by the surface area of the corresponding cohesive element.
 Cohesive Energy 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 opening in the loading history.
 This item is not available in cohesive contact.
 Cohesive Opening(s)
 o This provides the output for openings in the 3 modes (Mode I, II, and III – which are aligned to the local elemental system) and the basic system.
 This item is not available in cohesive contact.
 Cohesive Status (Status)
 Indicates element loading/unloading/fail state
 0: Loading
 1: Unloading/Reloading
 2: Fail
 This item is not available in cohesive contact.
 Cohesive Traction (Traction)
 This provides the output for traction in the local elemental system (Mode I, II, and III) and the basic system.
 This item is not available in cohesive contact.
 Delamination Growth Index (Scalar)
 Represents the ratio of the Energy of the cohesive element divided by energy that the element can withstand before failure (that is, the element reaches the end of the tractionopening curve).
 This result is nonzero only when the damage starts to occur for a particular element.
 This item is not available in cohesive contact.
Eroded cohesive elements are not shown in h3d output starting from the eroded analysis time of the element. Therefore, the above items are not available in eroded cohesive elements. See Cohesive Element Erosion for more information.
When contacts are used to model the cohesive zone, the output is interpreted in the following way. 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.’