**1. Introduction**

The requirements for technical products are rising. Fulfilling these becomes more and more di fficult using mono material components. A hybrid component combines di fferent materials and can be tailored to meet specific requirements that might not be possible with a single material. The joining process to produce hybrid components from two or more mono material parts is usually at the end of the production chain. If a forming process is included, joining is either during or after the forming process. Doing the joining process before forming allows to treat the joining zone during forming geometrically and thermomechanically. This provides more flexibility to tailor the hybrid component regarding its requirements; this process is called Tailored Forming (Figure 1).

**Figure 1.** Tailored forming process chain [1].

The joining zone is intended to be firmly bonded. In the case of aluminium and steel, intermetallic compounds often arise. Their existence shows that a firm bond is obtained during joining. Subsequent heat treatments cause thickening of intermetallic compounds to about 2 μm [2]. Though practical investigations indicate advantages of joining zones without or with a non-measurable intermetallic compound thickness during forming [3].

Thin layers are often modelled using cohesive zone elements. Standard formulations use traction separation laws for the constitutive relation which are not su fficient here. Forming of hybrid bulk metal components might include severe membrane deformations that are additional deformation modes not describable by separations. Internal Thickness Extrapolation (InTEx) [4] and Membrane Mode Enhanced Cohesive Zone Elements (MMECZE) [5] are two novel approaches that consider additional loading directions. Internal Thickness Extrapolation aims on layers that are thin but not completely flat, reconstructs a volumetric geometry from a flat element topology and applies bulk material models. Membrane Mode Enhanced Cohesive Zone Elements are focused on flat joining zones and combine Traction Separation Laws for the separation loading directions with a bulk material capturing the damage contribution from membrane mode deformation of the joining zone. Here, the latter approach is more convenient due to the flatness of the joining zone.

#### **2. Membrane Mode Enhanced Cohesive Zone Elements**

Membrane Mode Enhanced Cohesive Zone Elements are extensively described in [5]. Here, key points are summarised that are necessary for the implementation and application.

The Continuum Damage Mechanics framework from Lemaitre [6] is a concept that can describe damage in various bulk material models. Damage (0 ≤ *D* ≤ 1) reduces the e ffective surface area for load transmission. Stress is distributed on a smaller portion of the material; the elastic and the plastic response are a ffected. For Membrane Mode Enhanced Cohesive Zone Elements the framework is applied at an interface. Instead of tensorial stresses, traction vectors transmit the load through the material. Figure 2 shows the joining zone with partial defects (∂*SD*) that reduce the surface area for load transmission (*t*). The e ffective traction due to a reduction of the load transmitting surface area is

**Figure 2.** Damage of an interface [5]..

The change of damage is related to the rate *D* which depends on the combination of the plastic multiplier . λ describing the amount of deformation and the strain energy density release rate *Y*. The latter is the amount of elastic energy stored in the material that is released due to damage softening

$$
\dot{D} = \dot{\lambda} \frac{1}{1 - D} \left(\frac{\mathcal{Y}}{\mathcal{S}}\right)^{s} \tag{2}
$$

(1)

$$\text{with} \quad Y = \frac{1}{2} \frac{1}{(1-D)^2} \left(\frac{t\_t^2}{E} + \frac{t\_s^2}{G}\right) \tag{3}$$

where *s* and *S* are damage related material parameters, *tt* and *ts* are the tension and shearing part of the traction, and *E* and *G* are the Young's and shearing modulus.

Traction Separation Laws (TSLs) are defined on the cohesive surface (black dashed in Figure 3) that describes the midplane of the initially coincident faces of the upper and lower component

**Figure 3.** Separation of a joining zone; initial configuration (**a**) and current configuration (**b**) [5].

TSLs describe the traction over separation relation. They are often curves with a nonlinear load increase that levels off smoothly for a rising displacement until failure. This might be realistic for the application of modelling adhesives. Within Tailored Forming brittle failure of an infinitesimal thin joining zone has to be described. The theoretically infinitesimal stiffness with sharp load release after failure is modelled using a large but finite artificial stiffness k to circumvent numerical problems

$$t = k \,\delta \tag{5}$$

where *t* is the traction resulting from the separation δ (see green vectors in Figure 3). This traction separation relation induces a reaction force in the direction contrary to the separation.

The membrane mode deformation can be modelled, see Figure 3, as angular and length changes of the vectors *A*ξ and *<sup>A</sup>*η in the initial configuration to *<sup>a</sup>*ξ and *a*η in the current configuration. Describing these vectors in a two-dimensional local coordinate system that follows the cohesive surface (two dimensional quantities with ˆ in the following) allows to compute a two-dimensional deformation gradient

$$F = \begin{array}{c} \frac{\partial \ell}{\partial \bar{X}} = \frac{\partial \ell}{\partial \bar{\zeta}} \cdot \frac{\partial \bar{\zeta}}{\partial \bar{X}} \\ \text{with } \frac{\partial \ell}{\partial \bar{\zeta}} = \left( \begin{array}{cc} \bar{A}\_{\bar{\zeta}} & \bar{A}\_{\eta} \end{array} \right) \text{and } \frac{\partial \bar{\zeta}}{\partial \bar{X}} = \left( \begin{array}{cc} \mathfrak{a}\_{\bar{\zeta}} & \mathfrak{a}\_{\eta} \end{array} \right)^{-1} \end{array} \tag{6}$$

Based on the deformation gradient a variety of material models can be applied; here the Lemaitre model is used. The plastic multiplier as amount of (plastic) deformation can be computed (right in Figure 4). Together with tractions from the traction separation relation (left) a damage rate ( . *D*) results. Accumulation of the damage *D* (or certainly high tractions *tn* and *ts*) might then cause interruption of the cohesion.

**Figure 4.** Merging cohesive zone modes and membrane modes [5].

The Membrane Mode Enhanced Cohesive Zone Element concept does not necessitate this, but failure is updated explicitly here. This means the load transmission is cut off in the time step after failure occurred. Such procedure is numerically very advantageous as failure induced instabilities are inhibited and the nonlinearity of the computation is limited to the nonlinearity of the forming process itself.

#### **3. Cohesive Modelling in MSC Marc**

Marc is a nonlinear finite element solver from MSC Software (MSC Marc Mentat 2019 Feature Pack 1, Newport Beach, CA, US). The pre- and postprocessor Mentat, also from MSC Software, allows to setup simulations for Marc. Linear and quadratic Cohesive Zone Elements are provided by the software package [7]. Figure 5 shows the connectivity (a) and Gauss integration points (b) of element 188. The connectivity corresponds to a brick element, though the thickness in one direction is (initially) zero so the geometry complies with a quadrilateral element. As integration is only executed in the cohesive surface (dashed midplane, Figure 5a), a Gauss integration scheme contains only 2 × 2 = 4 points like a quadrilateral element. Wedge or triangle like element types and Newton Cotes integration schemes are also available. The degrees of freedom are displacements in three directions. In case of a thermomechanical analysis a weak coupling is realized with a staggered approach, i.e., a separate thermal element is used additionally to the mechanical element (e.g., the thermal element 222 corresponds to the mechanical element 188).

**Figure 5.** MSC Marc element 188: connectivity (**a**) and Gauss integration points (**b**) [7].

Modelling of cohesive zones is realised as a subsequent step of solid modelling. Existing nodes are doubled and the Cohesive Zone Element is placed. Cohesive elements can either immediately be inserted (Toolbox/General/Matching Boundaries) or when a certain load criterion is reached (Toolbox/Fracture mechanics/Delamination). Here, the direct insertion is beneficial to track membrane changes even if the load remains small. The inserted cohesive elements can use user subroutines to define the material behaviour.

Figure 6 shows the analysis path from pre- to postprocessing; element and material subroutines are drawn as sub steps of processing to discuss the implementation that interacts here. During processing, coordinates and displacements are handed over to the element subroutine. In the case of cohesive elements, separations are computed and the material subroutine is called. The material subroutine returns a tangent and stress. This is processed in the element subroutine to an element stiffness matrix and a residual.

**Figure 6.** Processing.

#### **4. Membrane Mode Enhanced Cohesive Zone Element Technology as a Material Subroutine**

A cohesive formulation based on separations does not need to know more about the deformation than the separations. Though Membrane Mode Enhanced Cohesive Zone Elements need information about the membrane deformation. Therefore, in [5] they are implemented as a user element. To comfortably use this new technology within MSC Marc, it has to be implemented as a material model. The missing information has to be gathered and the differing output has to be discussed as well.

#### *4.1. Gathering Information about the Membrane Deformation*

For the sake of simplicity, Figure 6 only shows the separation as most important information that is transferred from the element to the material. Some further data are the external element number m(1) and the integration point number nn. However, the required information about the membrane deformation is not transferred. Though Marc provides some functions that help to gather the required information about the membrane deformation [8]. The (internal) element nodes nodes and the number of nodes num can be determined using

call elnodes(ielint(m(1)),num,nodes)

where ielint(m(1)) converts the external to an internal element number. The number of nodes is used to identify the current element type. The internal node numbers have to be converted to external node numbers

```
nodeid = nodext(nodes(i))
as these can be used to obtain nodal values with
call nodvar(icod,nodeid,valno,nqncomp,nqdatatype)
```
where, e.g., icod=0 gives coordinates and icod=1 displacements.

The requested data is returned in valno; nqncomp and nqdatatype are the size and type of the returned data. With these the initial and current cohesive surface geometry can be constructed; e.g. the corner points for a linear quad like element in initial configuration are

xs(:,1) = (xl(:,1) + xl(:,5))/2 xs(:,2) = (xl(:,2) + xl(:,6))/2 xs(:,3) = (xl(:,3) + xl(:,7))/2 xs(:,4) = (xl(:,4) + xl(:,8))/2,

where xl are the node coordinates. For the current configuration the displacements have to be added. All this happens within a material subroutine that is called in one certain integration point. The location of integration point number nn for a linear quad with Gauss integration is

```
select case (nn)
  case (1)
    xi = − 1.d0/dsqrt(3.d0)
     eta = − 1.d0/dsqrt(3.d0)
  case (2)
    xi = + 1.d0/dsqrt(3.d0)
     eta = − 1.d0/dsqrt(3.d0)
  case (3)
    xi = − 1.d0/dsqrt(3.d0)
     eta = + 1.d0/dsqrt(3.d0)
  case (4)
    xi = + 1.d0/dsqrt(3.d0)
     eta = + 1.d0/dsqrt(3.d0)
end select .
With standard bilinear ansatz functions and derivativesN(1) = 0.25d0 * ( 1.d0 − xi) * ( 1.d0 − eta)
N(2) = 0.25d0 * ( 1.d0 + xi) * ( 1.d0 − eta)
N(3) = 0.25d0 * ( 1.d0 + xi) * ( 1.d0 + eta)
N(4) = 0.25d0 * ( 1.d0 − xi) * ( 1.d0 + eta)
Nxi(1) = 0.25d0 * (−1.d0 + eta)
Nxi(2) = 0.25d0 * ( 1.d0 − eta)
Nxi(3) = 0.25d0 * ( 1.d0 + eta)
Nxi(4) = 0.25d0 * (−1.d0 − eta)
Neta(1) = 0.25d0 * (−1.d0 + xi)
Neta(2) = 0.25d0 * (−1.d0 − xi)
Neta(3) = 0.25d0 * ( 1.d0 + xi)
Neta(4) = 0.25d0 * ( 1.d0 − xi)
```
the 3D tangential vectors in initial and current configuration can be determined. E.g., for the quadrilateral like element with the number of cohesive surface points num2 the tangential vectors in the initial configuration are

```
axi3Di = (/ 0.d0, 0.d0, 0.d0 /)
aet3Di = (/ 0.d0, 0.d0, 0.d0 /)
do i = 1, num2
  axi3Di(1) = axi3Di(1) + xs(1,i) * Nxi(i)
  axi3Di(2) = axi3Di(2) + xs(2,i) * Nxi(i)
  axi3Di(3) = axi3Di(3) + xs(3,i) * Nxi(i)
  aeta3Di(1) = aeta3Di(1) + xs(1,i) * Neta(i)
  aeta3Di(2) = aeta3Di(2) + xs(2,i) * Neta(i)
  aeta3Di(3) = aeta3Di(3) + xs(3,i) * Neta(i)
```