# Thermomechanical modeling of copper-ceramic bonded assemblies for power electronics applications

Théophile Hourdou

► To cite this version:

Théophile Hourdou. Thermomechanical modeling of copper-ceramic bonded assemblies for power electronics applications. Univ. Grenoble Alpes, CNRS, Grenoble INP, SIMaP, 38000 Grenoble, France. 2024. hal-04682364

HAL Id: hal-04682364

https://hal.science/hal-04682364v1

Submitted on 30 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Thermomechanical modeling of copper-ceramic bonded assemblies for power electronics applications

Théophile Hourdouª

ªUniv. Grenoble Alpes, CNRS, Grenoble INP, SIMaP, 38000 Grenoble, France

August 30, 2024

This work is supported by the French National Research Agency in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02)

**Abstract:** This study examines the thermo-mechanical behaviour and reliability of Direct Bonded Copper (DBC) component in power electronics. Using finite element simulations and experimental data, material properties and interface behaviour has been modelled. Thermal cycling simulation has been done to identify key failure modes, focusing on stress distribution and crack propagation at the copper-ceramic interface. The findings offer crucial insights for improving the design and durability of DBC components in demanding applications.

## 1 Introduction

In the field of power electronics, the demand for highly reliable and efficient thermal management solutions has significantly increased, driven by the rapid development of high-power density applications such as electric vehicles, renewable energy systems, and advanced computing technologies. Direct Bonded Copper (DBC) (Figure 1) components have emerged as a critical one in these applications, offering a combination of high thermal conductivity and electrical insulation in packaging devices.

<figure><img src="image_2.png" /><figcaption>Figure 1: Global (a) and sectional view (b) of a power electronics device

[Ben Kaabar, 2015]</figcaption></figure>

However, the industrial adoption of DBC components faces a critical challenge due to their susceptibility to failure under thermal cycling conditions. The mismatch in the coefficients of thermal expansion (CTE) between the copper and ceramic layers introduces significant thermo-mechanical stresses during operation, which can lead to the degradation of the material interfaces and ultimately result in device failure.
<figure><img src="image_3.png" /><figcaption>(a) crack propagation in Al₂O₃ ceramic</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) ceramic-copper interface decohesion</figcaption></figure>

<caption>Figure 2: Failures in a DBC component [Gaiser et al., 2015]</caption>

Figure 2 shows the two primary failure mechanisms of the DBC component. A common one is crack initiation and propagation in the underneath ceramic. However, interface decohesion, where the bond between the copper and ceramic layers weakens and separates, is the more prevalent failure mode.

## 1.1 Objectives

Understanding the behaviour of materials and their interfaces under these conditions is paramount. This study is focused on modelling approaches to predict the main failure mechanisms being ceramic-copper interface decohesion. This involves not only the precise mechanical characterization of the individual materials but also the development of fine numerical models that can simulate the complex interactions between the copper, ceramic, and their interfaces. The outcomes of this study are expected to contribute to the advancement of power electronics with more reliable DBC components.

## 1.2 Methodology

This study uses experimental data from [Ben Kaabar, 2015] to model the thermo-mechanical behaviour of DBC components and interface. Two methodologies has been used to analyze the behaviour and reliability of the DBC component.

First, a simulation of the DBC component cooling during assembly is conducted to quickly analyse the interactions between the materials under operational conditions. The study then focuses on evaluating interface failure with the coupled criterion, using the analysis of stress levels along the interface and the energy release rate associated with crack propagation, providing insights into the conditions under which cracks may initiate at the interface.

Secondly, a simulation of the DBC component during assembly and thermal cycling during operation is conducted, using the thermo-mechanical non-linear behaviours of the materials and the failure behaviour of the interface. The study then focuses on evaluating interface damage and failure.
## 2 DBC components behaviour modeling

For both methodologies, thermo-mechanical model behaviour parameters must be identified at first from experimental data. The global behaviour of bulk components is based on the total strain $\underline{\epsilon}$ additive decomposition:

$$ \underline{\epsilon} = \underline{\epsilon}_e + \underline{\epsilon}_p + \underline{\epsilon}_{th} \quad (1) $$

with subscript e the reversible (or elastic) strain part, p the irreversible (or plastic) strain part if appropriate and th the thermal strain part.

**Reversible strain part** is commonly modelled by isotropic linear elasticity (Hooke's law):

$$ \underline{\epsilon}_e = \frac{1 + \nu}{E} \underline{\sigma} - \frac{\nu}{E} I_1(\underline{\sigma}) \underline{I} \quad (2) $$

where E is the Young's modulus and ν the Poisson's ratio.

**Irreversible strain part** in metallic materials is commonly described by an equivalent stress $\sigma_{eq}$ determined from a yield criterion and various internal variables.

The equivalent stress is defined as follows:

$$ \sigma_{eq}(\sigma_{*}, R) = \max\{\sigma_{*} - R, 0\} \quad (3) $$

with $\sigma_{*}$ an effective stress determined from the yield criterion and R an isotropic hardening. In the literature, this effective stress is determined by the von Mises yield criterion:

$$ \sigma_{*} = J_{2} \left( \underline{\sigma} - \sum_{i=1}^{n} \underline{X}_{i} \right) \quad (4) $$

with $\underline{X}$ kinematic hardenings.

Several evolution functions may be used to describe both hardening, however in the present study non-linear ones are used since metallic behaviour exhibit severe non-linearities. The isotropic hardening may be defined as a function of the cumulated plastic strain $p$:

$$ R = R_{0} + Q(1 - \exp(-bp)) \quad (5) $$

with $R_0$ the initial yield stress, Q the hardening amplitude and b the saturation rate.

Kinematic hardenings are defined by their rate, depending on the plastic strain rate $\dot{\underline{\epsilon}}_p$:

$$ \dot{\underline{X}}_i = \frac{2}{3} C_i \dot{\underline{\epsilon}}_p - \gamma_i \underline{X}_i \dot{p} \quad (6) $$

with $C_i$ the hardening rate and $\gamma_i$ the saturation rate.

**Thermal strain part** is supposed isotropic and depends linearly on the temperature:

$$ \underline{\epsilon}_{th} = \alpha \Delta T \underline{I} \quad (7) $$

with $\alpha$ the coefficient of thermal expansion (CTE).

This coefficient was measured experimentally, values were used by [Ben Kaabar, 2015] and are kept as is in this study. Furthermore, an assumption is made about the temperature independence of all mechanical parameters i.e. all values are considered constant.
## 2.1 Ceramic thermo-mechanical modeling

Ceramic material exhibits a linear response under traction and bending as shown by previous studies [Ben Kaabar, 2015]. Experiments results have already been used to finely model its mechanical behaviour. Results and parameters are only presented here as a recall.

### 2.1.1 Experimental data

3-point bending were used to characterize the ceramic material. A parallelepiped-shape sample is placed on 2 lateral pins and under a third central one. Figure 3(a) shows the experimental setup. Typical results of ceramic bending, which exhibit a linear evolution of the load depending on the pins displacement up to failure, are shown in Figure 3(b).

<figure><img src="image_5.png" /><figcaption>Figure 3: Ceramic mechanical characterization</figcaption></figure>

### 2.1.2 Model description

Since experimental results are linear, the material mechanical behaviour may be modelled by a linear elastic model presented in the previous section, in equation (2). The complete thermo-mechanical model is defined with the following material parameters:

* elasticity : E and ν

* thermal : α

### 2.1.3 Parameters identification

[Ben Kaabar, 2015] study used beam theory to recover E from the following equation:

$$E = \frac{F_R L^3}{48 \delta I_z} \qquad (8)$$

with $F_R$ the load at failure, $L$ the length between lateral pins, $δ$ the deflection (here supposed equal to the central pin displacement) and $I_z$ the second moment of area of the beam's cross section depending on the section geometry.

Concerning Poisson's ratio, [Ben Kaabar, 2015] study evaluated its value by homogenization and comparison to the literature.
Final parameters values shown in Table 1 were found to modelled the ceramic material.

<table><caption>Table 1: Ceramic thermo-mechanical model parameters</caption><thead><tr><th>E [GPa]</th><th>ν [-]</th><th>α [e-6/°C]</th></tr></thead><tbody><tr><td>330</td><td>0.22</td><td>8</td></tr></tbody></table>

## 2.2 Copper thermo-mechanical modeling

Unlike ceramic in the previous section, copper exhibit a non-linear response under cyclic tensile tests as shown by previous studies [Ben Kaabar, 2015] and was modelled by an elasto-plastic mechanical model with a single cinematic non-linear hardening. In this study, the model has been re-identified to be more representative of experimental data<fnref n="1" />.

### 2.2.1 Experimental data

Cyclic tensile tests were done in [Ben Kaabar, 2015] study using plate dogbone samples. 10 tensions up to a load corresponding to an engineering stress of 120MPa was imposed. Load cell were used to acquire the load and digital image correlation were used to compute longitudinal and transverse extensions in the calibrated length of the sample.

True stress and true strain were computed from load and extensions. Figure 4 shows the resulting cyclic tensile curve for this metallic material. Left side corresponds to the total test, right side shows a zoom around 5% of strain corresponding to additional cycles.

<figure><img src="image_6.png" /><figcaption>Figure 4: Cyclic tensile test data</figcaption></figure>

During the first tension, the material exhibits a linear response up to ~50MPa followed by a short severely non-linear curve (up to ~65MPa), then followed by a relatively long non-linear curve up to 120MPa.

Following cycles show that the severe non-linear curve occurs at each tension. Moreover, this non-linearity tends to start at increasingly stress values cycle after cycle, leading to the closing of the hysteresis loops and the stabilization of the mechanical response.

<footnote n="1">¹It may be noted that the data has been recovered graphically from [Ben Kaabar, 2015] thesis since numerical datasets had been lost.</footnote>
### 2.2.2 Model description

Since experimental results are non-linear and exhibit reversible and irreversible strains, the material mechanical behaviour may be modelled by an elasto-plastic model presented in previous section.

Elastic part is modelled with Hooke's law in equation (2).

Plastic part is modelled with von Mises criterion and several hardening functions. Two kinematic hardenings are chosen since a short and severe non-linearity occurs at each cycle and a relatively long one is observed in the first one. Moreover, an isotropic hardening is chosen to represent the stress value increase at each cycle.

The complete thermo-mechanical model is defined with the following material parameters:

* elasticity: $E$ and $\nu$

* plasticity: $R_0$, $C_1$ and $\gamma_1$, $C_2$ and $\gamma_2$, $Q$ and $b$

* thermal: $\alpha$

### 2.2.3 Parameters identification

Identification of parameters were done directly from the cyclic tensile curve.

Initial elastic slope allow determination of $E$. The plastic strain is then determined by removing the elastic part: $\epsilon_p = \epsilon - \sigma/E$. From the chart $\sigma$ vs. $\epsilon_p$, initial yield stress and hardening parameters are roughly evaluated.

An optimization by inverse method is used afterwards to refine the relevance of these parameters. It consists of a minimization of the deviation between experimental and simulated data. Identification results are shown in Figure 5. An excellent agreement between experimental data and numerical results can be seen.

<figure><img src="image_7.png" /><figcaption>Figure 5: Copper cyclic tensile curve</figcaption></figure>

Poisson's ratio and CTE are determined from literature.
Final parameters values shown in Table 2 were found to modelled the copper material.

<table><caption>Table 2: Copper thermo-mechanical model parameters</caption><thead><tr><th>E [GPa]</th><th>ν [-]</th><th>α [-/°C]</th><th>R<sub>0</sub> [MPa]</th><th>C<sub>1</sub> [MPa]</th><th>γ<sub>1</sub> [-]</th><th>C<sub>2</sub> [MPa]</th><th>γ<sub>2</sub> [-]</th><th>Q [MPa]</th><th>b [-]</th></tr></thead><tbody><tr><td>127</td><td>0.33</td><td>16.5e-6</td><td>43</td><td>27000</td><td>1850</td><td>1120</td><td>15</td><td>30</td><td>12.8</td></tr></tbody></table>

## 2.3 Interface mechanical modeling

Last part of the BDC component is the interface between ceramic and copper layers. Interfaces must be modelled differently from bulk materials with an approach linked to fracture mechanics.

When the interface is subjected to relatively low stress values, neither crack occurs nor decohesion propagation. The interface simply represents a discontinuity between to homogeneous sections.

With higher stress levels, the interface begins to open and a degradation of its mechanical properties occurs, at the decohesion front, prior to crack opening. This behaviour is well modelled by the cohesive zone method with a common and simple bilinear one of [Dugdale, 1960] shown in Figure 6.

<figure><img src="image_8.png" /><figcaption>Figure 6: Cohesive zone method [Ben Kaabar, 2015]</figcaption></figure>

In Figure 6 left side, a crack inside a plate under tension is shown. The half-crack length is noted *a* from its symmetry axis. In this area, crack lips are stress free. Further, a stress field is however visible along a short length even if the crack seems to open in this region, meaning a separation behaviour with non-zero rigidity.

In Figure 6 right side, the proposed separation behaviour is shown with normal stress vector component $T_n$ vs. crack opening $\delta_n$. Two stress evolution phases occur during opening. First, the crack is allowed to slightly open reversibly up to a maximum stress value $\sigma_c$ with a rigidity $k = \sigma_c/\delta_0$. Secondly, damage initiation occurs, and increase, up to a maximum opening at failure $\delta_{cr}$. During this second phase, the material weakens and the stress decrease linearly. The energy used to open the crack $\Phi_{sep}$ is simply the integration of the two phases on the figure, also named $G_c$ as the energy release rate.
### 2.3.1 Experimental data

In order to characterized interface decohesion in DBC component between ceramic and copper layers, 4-point bending tests of these components were done in [Ben Kaabar, 2015] study. A parallelepiped sample with a central layer of ceramic and two layers of copper were bent between 2 inner bottom pins and 2 outer top pins. A notch is machined in the upper copper layer to ensure ceramic crack between inner pins. As seen on Figure 7, the displacement $\Delta L$ of the outer ones is imposed and the reaction load $F$ is measure on the inner ones.

<figure><img src="image_9.png" /><figcaption>Figure 7: DBC bending test</figcaption></figure>

The test is instrumented with an optical camera in order to observe the decohesion propagation during bending. Figure 8 shows the load vs. displacement with optical data at specific time during the test. Two phases are visible. A first one corresponding to the bending of the undamaged structure up to point A where the ceramic layer break inside the notch region. At this point, decohesion initiation occurs between ceramic layer and the underneath copper one. Phase B corresponds to the propagation of the decohesion in both direction. However for this specific test, it has been seen a longer decohesion at the right side leading to a significant di-symmetry.

<figure><img src="image_10.png" /><figcaption>Figure 8: DBC bending test results</figcaption></figure>

A processing of the optical data were done by [Ben Kaabar, 2015] on this right side to evaluate the decohesion length. However an additional analysis were done in the present study to compare left and right propagation rates. The results are shown in the following sections.
### 2.3.2 Model description

The [Dugdale, 1960] model is well defined for interface behaviour during decohesion between two different materials. Therefore the complete mecanical model is defined with the following interface parameters based on Figure 6:

* rigidity $k$

* maximum stress $\sigma_c$

* energy release rate at failure $G_c$

### 2.3.3 Parameters identification

An optimization by inverse method is used to identify these interface model parameters. It consists of a minimization of the deviation between experimental and simulated data. A simulation of the 4-point bending test was done by finite-element method using Simulia Abaqus software. The simulation input script is available in appendix A. The structure was meshed and a partial view, including the initial crack in the ceramic layer and mesh refinement along interface with underneath copper layer, can be seen in Figure 9.

<figure><img src="image_11.png" /><figcaption>Figure 9: Simulation mesh of the DBC bending test</figcaption></figure>

Identification results are shown in Figure 10 with a comparison of load over sample width vs. pin displacement with left and right decohesion length as well. An excellent agreement between experimental data and numerical results can be seen.

<figure><img src="image_12.png" /><figcaption>Figure 10: Interface identification results</figcaption></figure>

In order to model well the di-symmetry, left and right interfaces have been defined with different parameters values. Results can be seen in Table 3. The small difference in these values is enough to generate such a propagation di-symmetry.
<table><caption>Table 3: Interface mechanical model parameters</caption><thead><tr><th></th><th>k [N/mm]</th><th>σ<sub>c</sub> [MPa]</th><th>G<sub>c</sub> [N/mm]</th></tr></thead><tbody><tr><td>left</td><td>1e8</td><td>160</td><td>0.120</td></tr><tr><td>right</td><td>1e8</td><td>150</td><td>0.115</td></tr></tbody></table>

# 3 Simulation of axi-symmetric DBC cooling

After material modelling, the structure can be modelled and simulated as well. In this study, finite element method was used to simulate temperature cooling step during DBC assembly which may involve interface decohesion due to CTE mismatch between ceramic and copper layers.

## 3.1 Finite element model

A finite element model was developed with simplified geometry and material behaviours. In this model a local part of the DBC structure is considered using Simulia Abaqus software. The simulation input script is available in appendix B.

<figure><img src="image_13.png" /><figcaption>Figure 11: Geometry (a) and mesh (b) of the DBC structure</figcaption></figure>

Figure 11 shows the structure with symmetry assumption and the two used materials. Kinematic boundary conditions are imposed:

1. axi-symmetry along left dashed line - U1 = 0

2. vertical symmetry along bottom dotted line - U2 = 0

Temperature field is imposed homogeneous in the structure, decreasing from 1000°C to 0°C.

Copper and ceramic are modelled as linear elastic materials. Only Young's modulus, Poisson's ratio and CTE are defined in the model. This model is indeed a first evaluation of interface decohesion since a quick and simple determination of stress levels and energy release rate is made, as defined by fracture mechanics.

Stress levels at the end of cooling are directly evaluated at interface nodes. The stress
field is extrapolated at nodes and the stress vector is computed:

$$ \underline{T} = \underline{\sigma n} \qquad (9) $$

with $n$ the norm of the interface at the node. Since the interface is oriented along the X-axis, the stress vector is directly known from the stress tensor components.

However energy release rate needs additional computation, here by the node release method. The interface nodes are duplicated and linked with kinematic equations. Virtual propagation steps are simulated by deactivating these equations node after node, starting from the right corner toward the axi-symmetry axis. The global strain energy of the structure $W_e$ is computed and extracted for each propagation step. The energy release rate $G$ is then defined as:

$$ G(a) = \frac{W_e(a) - W_e(a = 0)}{S(a)} \qquad (10) $$

with $a$ the length of the crack newly created and $S$ the surface of the crack.

## 3.2 Results

Simulations of several geometries were done to evaluate their impact on $\underline{T}$ and $G$. At first, copper thickness, corresponding to 1/4 the ceramic thickness initially, was multiplied by a factor of x2, x4, x8 and x16. The results are shown in Figure 12.

<figure><img src="image_14.png" /><figcaption>Figure 12: Axi-symmetric simulation results with copper thickness variation</figcaption></figure>

From Figure 12(a), it may be detrimental to increase copper thickness since the stress levels increase. However from Figure 12(b), $G$ levels increase up to a thickness ratio of 1/1 and then decrease. Therefore, this configuration where identical ceramic and copper layer must not be used.

A second morphology evaluation was done. Copper layer can be machined to add or remove volume of material. Holes (also called dimples) are created by removing copper near the corner. Therefore 4 dimples are meshed and are shown in Figure 13.
<figure><img src="image_15.png" /><figcaption>Figure 13: Dimple geometries</figcaption></figure>

Stress levels and energy release rate results are shown in Figure 14. Creating a deeper hole leads to reduce stress levels, however they increase significantly after the hole. The most important result is shown in Figure 14(b) where this dimple highly reduces the levels of G globally.

<figure><img src="image_16.png" /><figcaption>Figure 14: Axi-symmetric simulation results with dimples</figcaption></figure>

Additional morphology variations may be roughly evaluated based on this study. It has shown its capability to evaluate this criticity using a simple finite element simulation of an axi-symmetric structure with elastic components.

# 4 Simulation of 3D DBC cooling

To evaluate the decohesion propagation in a finer way, the model has to be more developed. Here the non-linear behaviour of the copper is used and an actual interface is modelled with cohesive zone method in the finite element calculation.

## 4.1 Finite element model

The structure is here 3D as shown in Figure 15. Symmetry plane at the structure back left, back right and bottom are taken into account.
<figure><img src="image_17.png" /><figcaption>Figure 15: Geometry of the DBC structure</figcaption></figure>

Symmetry boundary conditions are imposed and temperature field is imposed homogeneous in the structure, decreasing from 1000°C to 0°C. 20 heat-cooling cycles are simulated afterwards from -50°C to 250°C. The simulation was created using Simulia Abaqus software. Its input script is available in appendix C.

## 4.2 Results

Using cohesive zone model, decohesion is supposed to initiate and propagate during thermal cycles.

Figure 16 shows equivalent stress field after cooling (a) and after 20 cycles (b). The field is displayed on the deformed configuration. Stress level significantly decrease after cycles but no re-localisation can be clearly seen.

<figure><img src="image_18.png" /><figcaption>Figure 16: 3D simulation equivalent stress</figcaption></figure>

The deformed configuration shows an opening of the interface at the corner of the structure. Moreover, a comparison of interface damage is shown in Figure 17 between these two steps. The propagation is clearly visible with a concentration at the corner but also along copper edge.
<figure><img src="image_19.png" /><figcaption>Figure 17: 3D simulation interface damage</figcaption></figure>

From these results, this simulation has shown its capability to evaluate the propagation
in a quantitative way using finite element method.

**5 Conclusion**

This study provides significant insights into the thermo-mechanical behaviour and reliability of Direct Bonded Copper (DBC) components, which are critical components in power electronics. Through comprehensive modelling and simulation techniques, including the analysis of material behaviours, interface integrity, and the impact of thermal cycling, this study identifies key factors that contribute to DBC failure. The finite element simulations, coupled with experimental data, offer a detailed understanding of stress distribution and failure mechanisms, particularly at the copper-ceramic interface. The findings underscore the importance of accurate modelling and advanced design strategies to enhance the durability of DBC substrates under demanding operational conditions. This work not only advances the current knowledge of DBC component performance but also provides a robust framework for future studies aimed at improving the reliability of power electronics components.
A Source code of 4-point bending simulation

import numpy as np

from abaqus import mdb as _mdb, session as _session

import part as _part

import mesh as _mesh

import section as _section

import material as _material

import step as _step

import interaction as _interaction

import visualization as _visualization


class geometry:
    sample_length = 55. # mm: sample total length
    crack_position = 0.5 # mm: position from symmetry plane (mid-length) of the initial crack in the ceramic layer
    crack_angle = 90. #degree
    decohesion_length = 0.02 # mm: length of initial decohesion length from crack end

class ceramic:
    thickness = 0.25 # mm
    young = 330000. #GPa
    poisson = 0.22 # -

class copper:
    thickness = 0.5 # mm
    notch = 2. # mm: notch length
    young = 127000. #GPa
    poisson = 0.33 # -
    R0 = 43.2 # MPa
    C1 = 27000. # MPa
    D1 = 1850. # -
    C2 = 1120. # MPa
    D2 = 14.8 # -
    Q = 30. # MPa
    b = 12.8 # -

class _interface:
    length = 1.5 # mm
    k = 1.e8 # N/mm
    Tmax = 150. # MPa
    Gc = 0.01 # N/mm
    interfaces = (_interface(), _interface())

class pins:
    diameter = 4. # mm
    outer_span = 50. # mm
    inner_span = 25. # mm
    inner_lever = 20. # mm

class mesh:
    global_size = 0.05 #0.1 # mm
    interface_size = 0.005 # 0.02 # mm

name = 'simulation'

model = None
```python
def create():
    '''
    Create the simulation model in memory.

    This function do not take any parameter but uses the module parameters
classes to generate the model.

    Example
    -------
    >>> import model
    >>> model.mesh_interface_size = 0.01 # Change some parameters
    >>> model.create()

    Warning
    -------
    The model is regenerated at each call using current module parameters
    '''
    # Delete previous generation and create a new one
    global _model_
    if _model_ is not None:
        _mdb.models.changeKey(fromName=_model_.name, toName='To-Delete')
    model = _mdb.Model(name=name)

    del _mdb.models['Model-1' if _model_ is None else 'To-Delete']

    # Create the base geometry
    sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)
    sketch rectangle(
        point1=(-geometry.sample_length/2., -copper.thickness),
        point2=(geometry.sample_length/2., ceramic.thickness+copper.thickness)
    )
    sample = model.Part(
        name='Sample',
        dimensionality=_part.TWO_D_PLANAR,
        type=_part.DEFORMABLE_BODY
    )
    sample.BaseShell(sketch=sketch)
    del model.sketches['sketch']

    # Create notch
    sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)
    sketch rectangle(
        point1=(-copper.notation/2., ceramic.thickness),
        point2=(copper.notation/2., ceramic.thickness+copper.thickness)
    )
    sample.Cut(sketch=sketch)
    del model.sketches['sketch']

    # Create partition
    alpha = _np.radians(geometry.crack_angle) if geometry.crack_angle != 90. else _np.nan
    crack_end = 0. if _np.isnan(alpha) else ceramic.thickness/_np.tan(alpha)
    regular_end = 0. if _np.isnan(alpha) else crack_end/ceramic.thickness*mesh_interface_size

    sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)
    sketch Line( # ceramic top left
```
point1=(-geometry.sample_length/2., ceramic.thickness),
point2=(-copper.notation/2., ceramic.thickness)
)

sketch LINE(
  # ceramic top right
  point1=(copper.notation/2., ceramic.thickness),
  point2=(geometry.sample_length/2., ceramic.thickness)
)

sketch LINE(
  # ceramic bottom left
  point1=(-geometry.sample_length/2., 0.),
  point2=(geometry.crack_position-interfaces[0].length, 0.)
)

sketch LINE(
  # interfaces
  point1=(geometry.crack_position-interfaces[0].length, 0.),
  point2=(geometry.crack_position+interfaces[1].length, 0.)
)

sketch LINE(
  # ceramic bottom right
  point1=(geometry.crack_position+interfaces[1].length, 0.),
  point2=(geometry.sample_length/2., 0.)
)

sketch LINE(
  # underneath ceramic crack
  point1=(geometry.crack_position, -mesh Köhne interface size),
  point2=(geometry.crack_position, 0.)
)

sketch LINE(
  # underneath delamination left
  point1=(geometry.crack_position-geometry.decohesion_length, -mesh Köhne interface size),
  point2=(geometry.crack_position-geometry.decohesion_length, 0.)
)

sketch LINE(
  # underneath delamination right
  point1=(geometry.crack_position+geometry.decohesion_length, -mesh Köhne interface size),
  point2=(geometry.crack_position+geometry.decohesion_length, 0.)
)

sketch LINE(
  # ceramic crack
  point1=(geometry.crack_position, 0.),
  point2=(geometry.crack_position+crack_end, ceramic.thickness)
)

sketch LINE(
  # delamination left
  point1=(geometry.crack_position-geometry.decohesion_length, 0.),
  point2=(geometry.crack_position-geometry.decohesion_length+crack_end, ceramic.thickness)
)

sketch LINE(
  # delamination right
  point1=(geometry.crack_position+geometry.decohesion_length, 0.),
  point2=(geometry.crack_position+geometry.decohesion_length+crack_end, ceramic.thickness)
)

sketch LINE(
  # left interface end
  point1=(geometry.crack_position-interfaces[0].length, -mesh Köhne interface size),
  point2=(geometry.crack_position-interfaces[0].length, mesh Köhne interface size),
)

sketch LINE(
  # right interface end
  point1=(geometry.crack_position+interfaces[1].length, -mesh Köhne interface size),
  point2=(geometry.crack_position+interfaces[1].length, mesh Köhne interface size),
)

sketch LINE(
  # regular bottom
  point1=(geometry.crack_position-interfaces[0].length, -mesh Köhne interface size),
  point2=(geometry.crack_position+interfaces[1].length, -mesh Köhne interface size)
)

sketch LINE(
  # regular top
  point1=(geometry.crack_position-interfaces[0].length, mesh Köhne interface size),
  point2=(geometry.crack_position+interfaces[1].length, mesh Köhne interface size)
)
```pascal
)

sketch Line( # regular left
  point1=(geometry.crack_position-interfaces[0].length-mesh.global_size*3., -copper.thickness),
  point2=(geometry.crack_position-interfaces[0].length-mesh.global_size*3., ceramic.thickness)
)

sketch Line( # regular right
  point1=(geometry.crack_position+interfaces[1].length+mesh.global_size*3., -copper.thickness),
  point2=(geometry.crack_position+interfaces[1].length+mesh.global_size*3., ceramic.thickness)
)

sample.PartitionFaceBySketch(faces=sample.faces, sketch=sketch)

del model.sketches['sketch']

# Create pins
alpha=_np radians(30.)
radius = pins.diameter/2.-mesh.global_size
sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)
sketch.ArcByCenterEnds(
  center=(0.0, 0.0),
  point1=( pins.diameter/2.*_np.cos(alpha), pins.diameter/2.*_np.sin(alpha)),
  point2=(-pins.diameter/2.*_np.cos(alpha), pins.diameter/2.*_np.sin(alpha)),
  direction=_part.COUNTERCLOCKWISE
)

sketch.ArcByCenterEnds(
  center=(0.0, 0.0),
  point1=( radius*_np.cos(alpha), radius*_np.sin(alpha)),
  point2=(-radius*_np.cos(alpha), radius*_np.sin(alpha)),
  direction=_part.COUNTERCLOCKWISE
)

sketch Line(
  point1=(-radius*_np.cos(alpha), radius*_np.sin(alpha)),
  point2=(-pins.diameter/2.*_np.cos(alpha), pins.diameter/2.*_np.sin(alpha))
)

sketch Line(
  point1=( radius*_np.cos(alpha), radius*_np.sin(alpha)),
  point2=( pins.diameter/2.*_np.cos(alpha), pins.diameter/2.*_np.sin(alpha))
)

pin = model.Part(name='Pin', dimensionality=_part.TWO_D OccURAR, type=_part.DEFORMABLE_BODY)

pin.BaseShell(sketch=sketch)

del model.sketches['sketch']

# Get crack and interface edges
crack = sample.edges.findAt(
  ((geometry.crack_position+crack_end/2., ceramic.thickness/2., 0.), ), # crack
  ((geometry.crack_position+regular_end/2, mesh_interface_size/2., 0.), ), # regular ligament
  ((geometry.crack_position-geometry.decohesion_length/2., 0., 0.), ), # left decohesion
  ((geometry.crack_position+geometry.decohesion_length/2., 0., 0.), ) # right decohesion
)

bonds = (
  sample.edges.findAt((geometry.crack_position-interfaces[0].length/2., 0., 0.), ), # left interface
  sample.edges.findAt((geometry.crack_position+interfaces[1].length/2., 0., 0.), ), # right interface
)

# Create the assembly
assembly = model=rootAssembly
instance = assembly Instance(name='Sample', part=sample, dependent=True)
left_inner pin = assembly Instance(name='Left-Inner-Pin', part=pin, dependent=True)
right_inner pin = assembly Instance(name='Right-Inner-Pin', part=pin, dependent=True)
```
left_outer pin = assembly Instance (name='Left-Outer-Pin', part=pin, dependent=True)
right_outer pin = assembly Instance (name='Right-Outer-Pin', part=pin, dependent=True)
assembly rotate (
    instanceList= ( 'Left-Outer-Pin', 'Right-Outer-Pin' ),
    axisPoint=(0.0, 0.0, 0.0),
    axisDirection=(0.0, 0.0, 1.0),
    angle=180.0
)
assembly translate (
    instanceList= ( 'Left-Inner-Pin' , ),
    vector= (-pins inner span/2. , -copper.thickness-pins.diameter/2. , 0.0)
)
assembly translate (
    instanceList= ( 'Right-Inner-Pin' , ),
    vector= ( pins inner span/2. , -copper.thickness-pins.diameter/2. , 0.0)
)
assembly translate (
    instanceList= ( 'Left-Outer-Pin' , ),
    vector= (-pins outer span/2. , ceramic.thickness+copper.thickness+pins.diameter/2. , 0.0)
)
assembly translate (
    instanceList= ( 'Right-Outer-Pin' , ),
    vector= ( pins outer span/2. , ceramic.thickness+copper.thickness+pins.diameter/2. , 0.0)
)
# Mesh the parts
elem_type = _mesh ElemType (elemCode=_mesh.CPE4I, elemLibrary=_mesh.STANDARD)
sample.set (regions=(sample.faces,), elemTypes=(elem_type,))
sample.setMeshControls (regions=sample.faces, elemShape=_mesh.QUAD)
# Set global size
sample.seedPart (size=mesh.global_size)
edges = sample_edges.findAt (
    ((-geometry.sample_length/4. , ceramic.thickness+copper.thickness, 0.), ),
    (( geometry.sample_length/4. , ceramic.thickness+copper.thickness, 0.), ),
    ((-geometry.sample_length/4. , ceramic.thickness, 0.), ),
    (( geometry.sample_length/4. , ceramic.thickness, 0.), ),
    (( 0. , -copper.thickness, 0.), ),
)
sample.seedEdgeBySize (edges=edges, size=mesh.global_size, constraint=_mesh.FINER)
# Set interface size
edges = sample_edges.findAt (
    (( geometry.crack_position-interfaces [0].length/2. , mesh_interface_size, 0.), ),
    (( geometry.crack_position-interfaces [0].length/2. , 0. , 0.), ),
    (( geometry.crack_position-interfaces [0].length/2. , -mesh_interface_size, 0.), ),
    (( geometry.crack_position+interfaces [1].length/2. , mesh_interface_size, 0.), ),
    (( geometry.crack_position+interfaces [1].length/2. , 0. , 0.), ),
    (( geometry.crack_position+interfaces [1].length/2. , -mesh_interface_size, 0.), ),
    (( geometry.crack_position+crack_end/2. -geometry.decohesion_length, ceramic.thickness/2. , 0.), ),
    (( geometry.crack_position+crack_end/2. , ceramic.thickness/2. , 0.), ),
    (( geometry.crack_position+crack_end/2. +geometry.decohesion_length, ceramic.thickness/2. , 0.), ),
    (( geometry.crack_position-mesh_interface_size/2. , -mesh_interface_size, 0.), ),
    (( geometry.crack_position-mesh_interface_size/2. , 0. , 0.), ),
    (( geometry.crack_position+regular_end-mesh_interface_size/2. , mesh_interface_size, 0.), ),
    (( geometry.crack_position+crack_end-mesh_interface_size/2. , ceramic.thickness, 0.), ),
    (( geometry.crack_position+mesh_interface_size/2. , -mesh_interface_size, 0.), ),
(( geometry.crack_position+mesh_interface_size/2., 0., 0.), ),
(( geometry.crack_position+regular_end+mesh_interface_size/2., mesh_interface_size, 0.), ),
(( geometry.crack_position+crack_end+mesh_interface_size/2., ceramic.thickness, 0.), ),
)
sample.seedEdgeBySize(edges=edges, size=mesh_interface_size, constraint=_mesh.FINER)

299

300 # Mesh the sample
301 sample(generateMesh())

302

303 # Mesh the pin
304 elem_type = _mesh.ElemType(elemCode=_mesh.CPE4R, elemLibrary=_mesh.STANDARD)
305 pin.setType(regions=(pin.faces,), elemTypes=(elem_type,))

306

307 pin.seedPart(size=mesh.global_size)
308 edges = pin.edges.findAt(
309 ((0., pins.diameter/2., 0.), ),
310 ((0., radius, 0.), ),
311 )

312 sample.seedEdgeBySize(edges=edges, size=mesh.global_size/2., constraint=_mesh.FINER)
313 pin(generateMesh())

314

315

316 # Create interface
317 etype = _mesh.ElemType(elemCode=_mesh.COH2D4, elemLibrary=_mesh.STANDARD, viscosity=0.001)

318

319 # Create crack: insert cohesive elements and delete them
320 nb_elements = len(samples)
321 region = _part Region(side1Edges=crack)
322 sample.insertElements(edges=region)

323

324 crack_set = sample.Set(name='Crack-Elements', elements=samples[nb_elements:] )
325 sample.setType(regions=crack_set, elemTypes=(etype,))

326

327 sample.deleteElement(sample elements[nb_elements:])
328 del sample sets['Crack-Elements']

329

330 # Insert cohesive elements in each interface
331 for side, interface in zip(’Left’, ’Right’), bonds):
332 nb_elements = len(sample elements)
333 region = _part .Region(side1Edges=_part.EdgeArray((interface,)) )
334 sample.insertElements(edges=region)
335
336 interface_set = sample. Set(name=side+’-Interface-Elements’, elements=sample .elements[nb_elements:] )
337 sample.setElementType(regions=interface_set, elemTypes=(etype,))

338

339

340 # Create sets
341 faces = sample faces. findAt(
342(-geometry.sample_length*0.45, ceramic.thickness/2., 0.), ),
343 (( geometry.crack_position-interfaces[0].length/2., ceramic.thickness/2., 0.), ),
344 (( geometry.crack_position+interfaces[1].length/2., ceramic.thickness/2., 0.), ),
345 (( geometry.crack_position+crack_end/2.+geometry.decohesion_length*2., ceramic.thickness/2., 0.), ),
346 (( geometry.crack_position+crack_end/2.-geometry.decohesion_length/2., ceramic.thickness/2., 0.), ),
347 (( geometry.crack_position+crack_end/2.+geometry.decohesion_length/2., ceramic.thickness/2., 0.), ),
348 (( geometry.sample_length*0.45, ceramic.thickness/2., 0.), ),
349 (( geometry.crack_position+regular_end/2.-geometry.decohesion_length/2., mesh_interface_size/2., 0.), ),
350 (( geometry.crack_position+regular_end/2.+geometry.decohesion_length/2., mesh_interface_size/2., 0.), ),
351 (( geometry.crack_position-interfaces[0].length/2., mesh_interface_size/2., 0.), ),
(( geometry.crack_position+interfaces[1].length/2., mesh_interface_size/2., 0.), )

) 352

) 353

ceramic_set = sample.Set(name='Ceramic-Set', faces=faces ) 354

355

356 faces = sample.faces.findAt(
  ((-geometry.sample_length*0.45, ceramic.thickness+copper.thickness/2., 0.), ), 357
  (( geometry.sample_length*0.45, ceramic.thickness+copper.thickness/2., 0.), ), 358
  ((-geometry.sample_length*0.45, -copper.thickness/2., 0.), ), 359
  (( geometry.sample_length*0.45, -copper.thickness/2., 0.), ), 360

361

362 (( geometry.crack_position-geometry.decohesion_length/2., -mesh_interface_size/2., 0.), ), 363
364 (( geometry.crack_position+geometry.decohesion_length/2., -mesh_interface_size/2., 0.), ), 365
366 (( geometry.crack_position+interfaces[0].length/2., -mesh_interface_size/2., 0.), ), 367
367 (( geometry.crack_position, -copper.thickness/2., 0.), ), 368
368 copper_set = sample.Set(name='Copper-Set', faces=faces ) 369

370 edge = sample.edges.findAt(
  ((-geometry.sample_length/4., ceramic.thickness+copper.thickness, 0.), ), 371
  ) 372
373 sample Surface(name='Top-Left-Surface', side1Edges=edge) 374

375 edge = sample.edges.findAt(
  (( geometry.sample_length/4., ceramic.thickness+copper.thickness, 0.), ), 376
  ) 377
378 sample Surface(name='Top-Right-Surface', side1Edges=edge) 379

380 edge = sample.edges.findAt(
  ((-geometry.sample_length*0.45, -copper.thickness, 0.), ), 381
  (( geometry.crack_position, -copper.thickness, 0.), ), 382
  (( geometry.sample_length*0.45, -copper.thickness, 0.), ), 383
  ) 384
385 sample Surface(name='Bottom-Surface', side1Edges=edge) 386

387 edge = pin.edges.findAt(( ( 0., pins.diameter/2., 0.), ), ) 388
388 pin Surface(name='Pin-Surface', side1Edges=edge) 389

390 pin Set(name='Pin', faces=pin faces) 391

392

393 # Assign materials
394 mat = model.Material(name='Pin') 395
396 mat.Elastic(table=((copper.young, copper.poisson), )) 396
397 model.HomogeneousSolidSection(name='Pin-Section', material='Pin', thickness=1) 397
398 pin SectionAssignment(region=pin sets['Pin'], sectionName='Pin-Section') 398

399 mat = model.Material(name='Ceramic') 400
401 mat.Elastic(table=((ceramic.young, ceramic.poisson), )) 401
402 model.HomogeneousSolidSection(name='Ceramic-Section', material='Ceramic', thickness=1) 402
403 sample SectionAssignment(region=ceramic_set, sectionName='Ceramic-Section') 403

404 mat = model.Material(name='Copper') 405
406 mat.Elastic(table=((copper.young, copper.poisson), )) 406
407 mat.Plastic(
  hardening=_material.COMBINED,
  dataType=_material.PARAMETERS,
  numBackstresses=2,
  table=((copper.R0, copper.C1, copper.D1, copper.C2, copper.D2), ) 407
408
409
410
```csharp
        )
    411
    412 mat.plastic.CyclicHardening parameters= True, table= ((copper.R0, copper.Q, copper.b), ) )
    413 model.HomogeneousSolidSection(name='Copper-Section', material='Copper', thickness=1)
    414 sample SectionAssignment(region=copper_set, sectionName='Copper-Section')
    415
    416 for side, interface in zip( 'Left', 'Right' ), interfaces:
    417     czm = model.Material(name=side+'-Interface')
    418     czm.Elastic( table= ((interface.k, interface.k, interface.k), ), type=_material.TRACTION)
    419     czm.QuadsDamageInitiation( table= ((interface.Tmax, interface.Tmax, interface.Tmax), ) )
    420     czm.quadsDamageInitiation.DamageEvolution( table= ((interface.Gc, ), ), type=_material.ENERGY)
    421
    422     model.CohesiveSection(
    423         name=side+'-Interface-Section',
    424         material=side+'-Interface',
    425         response=_section.TRACTION_SEPARATION
    426     )
    427     sample SectionAssignment(
    428         region=sample sets[side+'-Interface-Elements'],
    429         sectionName=side+'-Interface-Section'
    430     )
    431
    432
    433 # Create simulation step
    434 step = model StaticStep(
    435     name = 'Step',
    436     nlgeom = True,
    437     previous='Initial',
    438     initialInc=0.01,
    439     maxInc=0.01,
    440     maxNumInc=100000,
    441     minInc=1.e-16,
    442     timePeriod=4
    443 )
    444 step control.setValues(
    445     allowPropagation=False,
    446     resetDefaultValues=False,
    447     timeIncrementation=(4.0, 8.0, 9.0, 16.0, 10.0, 4.0, 12.0, 50.0, 6.0, 3.0, 50.0)
    448 )
    449 model.fieldOutputRequests[ "F-Output-1" ].setValues(
    450     variables=( 'U', 'RF', 'E', 'PE', 'PEEQ', 'S', 'SDEG', 'STATUS' ),
    451     timeInterval=0.05
    452 )
    453
    454
    455 # Create inner pins rigid body
    456 rpid = assembly REFERENCEPOINT( point=(0.0, -copper.thickness-pins.diameter/2. -pins inner lever, 0.0) ).id
    457 assembly. Set( referencePoints=(assembly referencePoints[ rpid], ), name='RefPoint')
    458
    459 edges = assembly. instances[ 'Left-Inner-Pin' ].edges + assembly. instances[ 'Right-Inner-Pin' ].edges
    460 model.R
    461     name='Constraint-1',
    462     refPointRegion= _part. Region( referencePoints= (assembly. referencePoints[ rpid], ) ),
    463     pinRegion= _part. Region( edges=edges)
    464 )
    465
    466
    467 # Create Boundary Conditions
    468 model.DisplacementBC(
    469     name='Inner-Pins',
```
createStepName='Initial',
region=assembly sets['RefPoint'],
u1=0.0,
u2=0.0
)
model.DisplacementBC(
name='Left-Outer-Pin',
createStepName='Step',
region=left_outer pin sets['Pin'],
u1=0.0,
u2=-4.0
)
model.DisplacementBC(
name='Right-Outer-Pin',
createStepName='Step',
region=right_outer pin sets['Pin'],
u1=0.0,
u2=-4.0
)

# Create frictionless contact
model.ContactProperty('Frictionless-Contact')
model.interactionProperties['Frictionless-Contact'] .NormalBehavior(
constraintEnforcementMethod=_interaction.AUGMENTED_LAGRANGE
)

# Create standard contact
model.ContactProperty('Contact')
model.interactionProperties['Contact'] .NormalBehavior(
constraintEnforcementMethod=_interaction.AUGMENTED_LAGRANGE
)
model.interactionProperties['Contact'] .TangentialBehavior(
formulation=_interaction.PENALTY,
table=((0.35, ), ),
fraction=0.005
)

model SurfaceToSurfaceContactStd(
name='Left-Inner-Contact',
createStepName='Initial',
master=instance surfaces['Bottom-Surface'],
slave=left_inner pin surfaces['Pin-Surface'],
sliding=_interaction SMALL, #FINITE,
enforcement=_interaction .NODE_TOok surface,
interactionProperty='Frictionless-Contact',
)
model SurfaceToSurfaceContactStd(
name='Right-Inner-Contact',
createStepName='Initial',
master=instance surfaces['Bottom-Surface'],
slave=right_inner pin surfaces['Pin-Surface'],
sliding=_interaction .SMALL, #FINITE,
enforcement=_interaction .NODE_TOok surface,
interactionProperty='Frictionless-Contact',
)
model SurfaceToSurfaceContactStd(
name='Right-Outer-Contact',
createStepName='Initial',
529 master=instance.surfaces['Top-Right-Surface'],
530 slave=right_outer.pin.surfaces['Pin-Surface'],
531 sliding=_interaction SMALL, #FINITE,
532 enforcement=_interaction. NODE_TOok surface,
533 interactionProperty='Frictionless-Contact',
534 )
535
536 model. SurfaceToSurfaceContactStd(
537 name='Left-Outer-Contact',
538 createStepName='Initial',
539 master=instance.surfaces['Top-Left-Surface'],
540 slave=left_outer.pin.surfaces['Pin-Surface'],
541 sliding=_interaction. SMALL, #FINITE,
542 enforcement=_interaction. NODE_TOok surface,
543 interactionProperty='Contact',
544 )
545
546
547 # Store the model globaly
548 _model_ = model
549
550
551
552 def update(side, Tmax, Gc):
553 ''''
554 Update simulation model interface parameters.
555
556 Parameters
557 -----------
558 side: str ('Left', 'Right')
559 Side to update
560
561 Tmax: float
562 New value of Tmax
563
564 Gc: float
565 New value of Gc
566
567 Example
568 -------
569 >>> import model
570 >>> model.create()
571 >>> model.update('left', 150, 0.01)
572 >>> model.write_input()
573
574 Raises
575 -------
576 RuntimeError if the model has not been generated
577 ''''
578 if _model_ is None:
579 raise RuntimeError('Model not generated')
580
581 interface = side+'-Interface'
582 interface = _model_.materials[interface]
583 interface.quadsDamageInitiation.setValues(table=((Tmax, Tmax, Tmax), ))
584 interface.quadsDamageInitiation.damageEvolution.setValues(table=((Gc, ), ))
585
586
587
```python
def write_input():
    '''
    Write the simulation inp file from the created model

592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
def post_process(odbname=None):
    '''
    Process simulation output data to compute decohesion length for each interface

615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
```

**Example**

```python
>>> from subprocess import run
>>> import model
>>> model.create()
>>> model.write_input()
>>> run('abaqus job=simulation interactive', shell=True)
>>> model.post_process()
```

**Example**

```python
>>> import model
>>> model.post_process('previous_run.odb')
```

**Raises**

```
RuntimeError if the model has not been generated and no output file is specified
'''

if __model__ is None and odbname is None:
    raise RuntimeError('Model not generated')

# Get a dedicated viewport
viewport = _session.viewports.values()[0]
viewport.makeCurrent()

# Open odb file
if odbname is None:
```
odbname = _model_.name
odb = _session.openOdb(odbname+'.'+odb', Premake = True, readInternalSets = True)
viewport.setValues (displayedObject = odb)

# Get load on inner pins reference point
data = _session .xyDataListFromField(
odb = odb,
outputPosition = _visualization .NODAL,
variable = (( 'RF' , _visualization .NODAL, (( _visualization .COMPONENT, 'RF2') ), )) ,
nodeSets = ( "REFPOINT" , )
)
results = _np .array(data[0])
with open( odbname+ '/load .sim', 'w') as fout:
fin = write( '#displacement [mm] load [N/mm] 'n')
 np .savetxt (fin, results)

# Get left decohesion length
instance = odb .rootAssembly .instances [ 'SAMPLE']
elements = instance.elementSets [ 'LEFT-INTERFACE-ELEMENTS'] .elements

nids = []
for element in elements:
nids. extend (element .connectivity[:2])
nids = _np .array(nids, type =int )-1
coordinates = _np .array([instance .nodes[nid] .coordinates[0] for nid in nids])
ranks = _np .argsort (coordinates)
coordinates = _np .unique (coordinates [ranks])

# - Get data at integration points
data = _session .xyDataListFromField(
odb = odb,
outputPosition = _visualization .INTEGRATION ,
variable = (( 'SDEG' , _visualization .INTEGRATION), ),
elementSets = ( "SAMPLE .LEFT-INTERFACE-ELEMENTS" , )
)
data = _np .array(data)[:, : ,1] [ranks]

# - Convert to single node position
crack = _np .empty(( len (coordinates), data .shape [1] ))
crack[0] = data [0]
crack[1:1] = (data [1:-2:2]+data [2:-1:2])/2.
crack[-1] = data [-1]

# - Get crack front position
for i, data in enumerate(crack.T):
for xa, xb, ya, yb in zip (coordinates [:-1], coordinates [1:], data [:-1], data [1:]):
if yb>0.999: break
if yb<0.1: results[i,1] =geometry .crack_position
elif yb>0.99999: results[i,1] =xb
else: results[i,1] =xa + (0.999-ya) * (xb-xa) / (yb-ya)

# - Compute crack length
results [1:1] = geometry .crack_position -results [1:1]
with open( odbname+ '/crack_left .sim', 'w') as fout:
fin .write( '#displacement [mm] length [mm] 'n')
 np .savetxt (fin, results)
```python
# Get right decohesion length
instance = odb.rootAssembly_instances['SAMPLE']
elements = instance.elementSets['RIGHT-INTERFACE-ELEMENTS'].elements

nids = []
for element in elements:
    nids.extend element.connectivity[:2])
    nids = _np.array(nids, dtype=int)-1
    coordinates = _np.array([instance.nodes[nid].coordinates[0] for nid in nids])
    ranks = _np argsort coordinates
    coordinates = _np.unique (coordinates [ranks]) [::: -1]

# - Get data at integration points
data = _session .xyDataListFromField(
    odb=odb,
    outputPosition= visualization .INTEGRATIONPOINT,
    variable=( 'SDEG', _visualization .INTEGRATION self),
    elementSets=("SAMPLE .RIGHT-INTERFACE-ELEMENTS", )
)
data = _np .array(data)[:,1] [ranks] [::: -1]

# - Convert to single node position
crack = _np .empty((len (coordinates), data .shape[1]))
crack[0] = data[0]
crack[1:-1] = (data[1:-2:2]+data[2:-1:2])/2.
crack[-1] = data[-1]

# - Get crack front position
for i, data in enumerate (crack.T):
    for xa, xb, ya, yb in zip (coordinates[:-1], coordinates[1:], data[:-1], data[1:]):
        if yb>0.999: break
        if yb<0.1: results[i,1] = geometry .crack_position
        elif yb>0.99999: results[i,1] = xb
        else: results[i,1] = xa + (0.999-ya) * (xb-xa) / (yb-ya)

# - Compute crack length
results [:,1] = geometry .crack_position
with open( odbname+ '/crack_right .sim ','w') as fout:
    fout.write([ mm] length[mm] \n')
    _np .savetxt (fout, results)
odb.close()
```
B Source code of axi-symmetric DBC simulation

1 import numpy as np
2
3 from abaqus import mdb as _mdb, session as _session
4 import part as _part
5 import mesh as _mesh
6 import section as _section
7 import material as _material
8 import step as _step
9 import load as _load
10 import visualization as _visualization
11
12
13
14 class ceramic:
15     length = 0.750 # mm
16     thickness = 0.317 # mm
17     young = 330000. #GPa
18     poisson = 0.22 # -
19     cte = 8.0e-6 # /dC
20
21 class copper:
22     length = 0.650 # mm
23     thickness = 0.127 # mm
24     young = 127000. # MPa
25     poisson = 0.33 # -
26     cte = 16.5e-6 # /dC
27
28 class mesh:
29     global_size = 0.05 #0.1 # mm
30     interface_size = 0.005 # 0.02 # mm
31
32 class loads:
33     class cooling:
34         initial = 1000.0 # dC
35         final = 0.0 # dC
36         duration = 1.0 # s
37
38 name = 'simulation'
39 _model_ = None
40
41
42
43 def create():
44     '''
45     Create the simulation model in memory.
46     '''
47     This function do not take any parameter but uses the module parameters
48     classes to generate the model.
49     '''
50 Example
51     >>> import model
52     >>> model.mesh.size = 0.01 # Change some parameters
53     >>> model.create()
54
55 Warning
56
```python
# Delete previous generation and create a new one

global _model_

if _model_ is not None:
    _mdb.models.changeKey(fromName=*_model_*.name, toName='To-Delete')

model = _mdb.Model(name=name)

del _mdb.models['Model-1' if _model_ is None else 'To-Delete']

# Create geometry

sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)

# Add a construction line for axisymmetry

sketch.ConstructionLine(point1=(0.0, -100.0), point2=(0.0, 100.0))

# Define contour lines

sketch LINE(
    point1=(0.0, 0.0),
    point2=(ceramic.length, 0.0)
)

sketch LINE(
    point1=(ceramic.length, 0.0),
    point2=(ceramic.length, ceramic.thickness)
)

sketch LINE(
    point1=(ceramic.length, ceramic.thickness),
    point2=(copper.length, ceramic.thickness)
)

sketch LINE(
    point1=(copper.length, ceramic.thickness),
    point2=(copper.length, ceramic.thickness+copper.thickness)
)

sketch LINE(
    point1=(copper.length, ceramic.thickness+copper.thickness),
    point2=(0.0, ceramic.thickness+copper.thickness)
)

sketch LINE(
    point1=(0.0, ceramic.thickness+copper.thickness),
    point2=(0.0, 0.0)
)

# Create the part

part = model.Part(
    name='DBC',
    dimensionality=_partAXISYMMETRIC,
    type=_part.DEFORMABLE_BODY
)

part.BaseShell(sketch=sketch)

del model.sketches['sketch']

# Create interface

sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)

sketch LINE(
    point1=(0.0, ceramic.thickness),
    point2=(copper.length, ceramic.thickness)
)

part.PartitionFaceBySketch(faces=part.faces, sketch=sketch)
```
116 del model.sketches['sketch']
117
118 # Create regular domain in copper
119 size = meshefa.nf.nf_size
120 sketch = model .ConstrainedSketch (name ='sketch', sheetSize =1.0)
121 sketch .Line(
122 point1 =(0.0, ceramic.thickness+size),
123 point2 =(copper.length, ceramic.thickness+size)
124 )
125 part .PartitionFaceBySketch(
126 faces =part .faces .findAt( (0.0, ceramic.thickness+copper.thickness, 0.0) ),
127 sketch =sketch
128 )
129 del model.sketches['sketch']
130
131 # Create regular domains in ceramic
132 sketch =model .ConstrainedSketch (name ='sketch', sheetSize =1.0)
133 sketch .Line(
134 point1 =(0.0, ceramic.thickness-size),
135 point2 =(copper.length, ceramic.thickness-size)
136 )
137 sketch .Line(
138 point1 =(copper.length, ceramic.thickness-size),
139 point2 =(copper.length, ceramic.thickness)
140 )
141 sketch .Line(
142 point1 =(copper.length, ceramic.thickness-size),
143 point2 =(copper.length+size, ceramic.thickness-size)
144 )
145 sketch .Line(
146 point1 =(copper.length+size, ceramic.thickness-size),
147 point2 =(copper.length+size, ceramic.thickness)
148 )
149 part .PartitionFaceBySketch(faces =part .faces .findAt( (0.0, 0.0, 0.0) ), sketch =sketch)
150 del model.sketches['sketch']
151
152 # Create reference points for Generalized Plane Strain
153 part .ReferencePoint(point =part .vertices .findAt( (0.0, 0.0, 0.0) ))
154
155 # Create assembly instance
156 assembly =model .rootAssembly
157 instance = assembly .Instance (name ='DBC', part =part, dependent =True)
158
159
160 # Setup the mesh
161 part .setMeshControls (regions =part .faces, elemShape = _mesh .QUAD)
162 etype = _mesh .ElemType (elemCode = _mesh .CAX4, elemLibrary = _mesh .STANDARD)
163 part .setElementType(regions =(part .faces,), elemTypes =(etype,))
164
165 # Set element sizes
166 part .seedPart (size =mesh .global_size, deviationFactor =0.01, minSizeFactor =0.1)
167
168 edge = part .edges .findAt( (copper.length/2, ceramic.thickness+copper.thickness, 0.0) )
169 part .seedEdgeBySize(edges =(edge,), size =mesh .global_size, constraint = _mesh .FINER)
170
171 edge = part .edges .findAt( (ceramic.length/2, 0.0, 0.0) )
172 part .seedEdgeBySize(edges =(edge,), size =mesh .global_size, constraint = _mesh .FINER)
173
174 edge = part .edges .findAt( (copper.length, ceramic.thickness+copper.thickness/2, 0.0) )
```python
175 part.seedEdgeByBias(
176     end1Edges=(edge,),
177     minSize=mesh_interface_size,
178     maxSize=mesh_global_size/2,
179     biasMethod=_mesh.SINGLE,
180     constraint=_mesh.FINER,
181 )
182
183 edge = part.edges.findAt((copper.length+ceramic.length)/2, ceramic.thickness, 0.0)
184 part.seedEdgeByBias(
185     end2Edges=(edge,),
186     minSize=mesh_interface_size,
187     maxSize=mesh_global_size/2,
188     biasMethod=_mesh.SINGLE,
189     constraint=_mesh.FINER,
190 )
191
192 interface = part.edges.findAt((copper.length/2, ceramic.thickness, 0.0))
193 part.seedEdgeBySize(edges=(interface,), size=mesh_interface_size, constraint=_mesh.FINER)
194
195 size = mesh_interface_size
196 edge = part.edges.findAt((copper.length/2, ceramic.thickness+size, 0.0))
197 part.seedEdgeBySize(edges=(edge,), size=mesh_interface_size, constraint=_mesh.FINER)
198
199 edge = part.edges.findAt((copper.length/2, ceramic.thickness-size, 0.0))
200 part.seedEdgeBySize(edges=(edge,), size=mesh_interface_size, constraint=_mesh.FINER)
201
202 # Generate mesh
203 part(generateMesh())
204 if part.getUnmeshedRegions():
205     raise RuntimeError('Mesh generation failed')
206
207
208 # Create bulk element sets
209 elements = part faces.findAt((0,0,0)).getElements()
210 elements += part faces.findAt((copper.length/2.0, ceramic.thickness-size/2.0, 0.0)).getElements()
211 elements += part faces.findAt((copper.length+size/2.0, ceramic.thickness-size/2.0, 0.0)).getElements()
212 part.Set(name='Ceramic-Elements', elements=elements)
213
214 elements = part faces.findAt((copper.length/2.0, ceramic.thickness+size/2.0, 0.0)).getElements()
215 elements += part faces.findAt((0, ceramic.thickness+copper.thickness, 0)).getElements()
216 part.Set(name='Copper-Elements', elements=elements)
217
218 # Generate interface elements
219 nb_elements = len(part)
220
221 region = _part Region(side1Edges=_part.EdgeArray((interface,)))
222 part.insertElements(edges=region)
223 interface_set = part.Set(name='Interface-Elements', elements=part[nb_elements:] )
224
225 etype = _mesh.ElemType(elemCode=_mesh.COHAX4, elemLibrary=_mesh.STANDARD, viscosity=0.001)
226 part.setType(regions=interface_set, etype=(etype,))
227
228 # Get nodes at interface linked to cohesive zone elements
229 nodes1, nodes2 = set(), set()
230 for element in interface_set:
231     nodes1.update element.connectivity[:2])
232     nodes2.update element.connectivity[2:])
233
```
```python
234 # Remove interface elements
235 part.deleteElement(interface_set)
236 del part.sets['Interface-Elements']
237
238 # Sort interface nodes by X position from end to symmetry plane
239 nodes1 = [part.nodes[nid] for nid in nodes1]
240 nodes2 = [part.nodes[nid] for nid in nodes2]
241 nid1 = _np.argsort([node coorinates[0] for node in nodes1])[:-1]
242 nid2 = _np.argsort([node coordinates[0] for node in nodes2])[:-1]
243 nodes1 = [nodes1[i] for i in nid1]
244 nodes2 = [nodes2[i] for i in nid2]
245
246 # Create node sets
247 part.Set(name='Nodes', elements=part.nodes)
248
249 nodes = [node for node in part.nodes if node coorinates[0] == 0]
250 part.Set(name="X0", nodes=_mesh.MeshNodeArray nodes])
251
252 nodes = _mesh.MeshNodeArray([node for node in part.nodes if node coorinates[1] == 0])
253 part.Set(name="Y0", nodes=nodes)
254
255 # Create interface node sets
256 part.Set(name='Interface-Nodes_Side1', nodes=_mesh.MeshNodeArray nodes1))
257 part.Set(name='Interface-Nodes_Side2', nodes=_mesh.MeshNodeArray nodes2))
258
259 # Create node sets for all release steps
260 for i, (node1, node2) in enumerate(zip nodes1, nodes2), start=1):
261     part.Set(name='Interface-Nodes_Side1'+str(i), nodes=_mesh.MeshNodeArray([node1]))
262     part.Set(name='Interface-Nodes_Side2'+str(i), nodes=_mesh.MeshNodeArray([node2]))
263
264
265 # Assign copper properties
266 mat = model.Material(name='Copper')
267 mat.Elastic(table=((copper.young, copper.poisson), ))
268 mat.Expansion(table=((copper.cte, ), ))
269
270 model.HomogeneousSolidSection(name='Copper-Section', material='Copper', thickness=copper.length)
271 part SectionAssignment(
272     region=part sets['Copper-Elements'],
273     sectionName='Copper-Section',
274 )
275
276 # Assign ceramic properties
277 mat = model.Material(name='Ceramic')
278 mat.Elastic(table=((ceramic.young, ceramic.poisson), ))
279 mat.Expansion(table=((ceramic.cte, ), ))
280
281 model.HomogeneousSolidSection(name='Ceramic-Section', material='Ceramic', thickness=ceramic.length)
282 part SectionAssignment(
283     region=part sets['Ceramic-Elements'],
284     sectionName='Ceramic-Section',
285 )
286
287
288 # Create steps
289 step = modelSlowSlowStep(
290     name='Cooling',
291     nlgeom=True,
302     previous='Initial',
```
293 initialInc=loads.cooling Duration/20.,
294 maxInc=loads.cooling duration/20.,
295 maxNumInc=int(1e6),
296 minInc=loads.cooling duration/1e30,
297 timePeriod=loads.cooling duration
298 )
299 step control.setValues(
300 allowPropagation=False,
301 resetDefaultValues=False,
302 timeIncrementation=(4.0, 8.0, 9.0, 16.0, 10.0, 4.0, 12.0, 50.0, 6.0, 3.0, 50.0)
303 )
304
305 model. fieldOutputRequests [ "F-Output-1" ]. setValues(
306 variables = ( 'U', 'E', 'PE', 'PEEQ', 'S', 'SDEG', 'STATUS' ),
307 timeInterval = loads.cooling duration/20.0,
308 )
309
310 num_steps = len (model parts [ 'DBC' ]. sets [ 'Interface-Nodes_Side1' ]. nodes )-1
311 previous = 'Cooling'
312 for i in range (1, num_steps+1):
313 step = 'Release'+str (i)
314
315 model. StaticStep(
316 name = step,
317 nlgeom = True,
318 previous = previous,
319 initialInc = 0.2,
320 maxInc = 0.2,
321 maxNumInc = int (1e6),
322 minInc = 1e-30,
323 timePeriod = 1.,
324 )
325 previous = step
326 model. fieldOutputRequests [ 'F-Output-1' ]. setValuesInStep (stepName ='Release1', timeInterval =1.0)
327 model. historyOutputRequests [ 'H-Output-1' ]. setValuesInStep (stepName ='Release1', timeInterval =1.0)
328
329
330 # Apply symmetry along X axis
331 model. DisplacementBC(
332 name ='X-symmetry',
333 createStepName = "Initial",
334 region = instance. sets [ 'X0' ],
335 u1 = _load. SET,
336 )
337
338 # Apply symmetry along Y axis
339 model. DisplacementBC(
340 name ='Y-symmetry',
341 createStepName = "Initial",
342 region = instance. sets [ 'Y0' ],
343 u2 = _load. SET,
344 )
345
346 # Apply temperature variation
347 model.Temperature(
348 name ='Cooling',
349 createStepName = "Initial",
350 region = instance. sets [ 'Nodes' ],
351 magnitudes =(loads.cooling initial, )
411 job.writeInput(consistencyChecking=False)
412 del _mdb.jobs[_model_.name]
413
414
415
416 def post_process(odbname=None):
417     '''
418     Process simulation output data to compute strain energy and stress vector along interface
419     Parameters
420         odbname: str, optional
421             When the model is not created, simulation output filename without extension
422     Example
423     '''
424     >>> from subprocess import run
425     >>> import model
426     >>> model.create()
427     >>> model.write_input()
428     >>> run('abaqus job=simulation interactive', shell=True)
429     >>> model.post_process()
430
431
432
433
434 Example
435     '''
436     >>> import model
437     >>> model.post_process('previous_run.odb')
438
439
440     '''
441     RuntimeError if the model has not been generated and no output file is specified
442     '''
443     if __model__ is None and odbname is None:
444         raise RuntimeError('Model not generated')
445
446     # Get a dedicated viewport
447     viewport = _session.viewports.values()[0]
448     viewport.makeCurrent()
449
450     # Open odb file
451     if odbname is None:
452         odbname = __model__.name
453         odb = _session.openOdb(odbname+'.'+odb', readonly=True, readInternalSets=True)
454     viewport.setValues displayObject=odb
455
456     # Remove field average to get contours by element
457     viewport.odbDisplay.basicConfig.VALUES(averageElementOutput=False)
458
459     # Use undeformed shape
460     viewport.odbDisplay.commonOptions.setValues(deformationScaling=_part.UNIFORM, uniformScaleFactor=0)
461
462     instance = odb.rootAssembly_instances['DBC']
463
464     # Get interface nodes
465     coordinates = [node cooridates[0] for node in instance.nodeSets['INTERFACE-NODES SIDE1'] .nodes]
466     ranks = _np.argsort( coordinates )[:-1]
467     crack_advance = _np.array(coordinates)[ranks]
468     interface_length = crack_advance[0]
469
crack_advance = interface_length-crack_advance

# Compute global system energy
energy = _session.XYDataFromHistory(
    name='Global strain energy',
    odb=odb,
    outputVariableName='Strain energy: ALLSE for Whole Model',
    steps=odb steps keys()
).data

energy = _np.array(energy)
energy = energy[energy[: ,0]>=1. ] # Remove cooling step
energy = energy[ _np.unique(energy[: ,0], return_index=True) [1] ] # Remove time duplicates
energy[: ,0] = crack_advance # Change from time to crack advance

# Save results to file
with open('global_energy.dat', 'w') as fout:
    fout.write('#crack_advance[mm] We[mJ]
')
 np savetxt (fout, energy)

# Clean data
del_session .xyDataObjects['Global strain energy']

# Compute stress vector components along interface
step = 'Cooling'
frame = odb .steps[step].frames[-1]
_session .odbData[odb.name].setValues(activeFrames=(( 'Cooling', (-1, )),))

# Get stress field for profiles
field = frame .fieldOutputs['S']

# Save nodal data in copper side
values = field.getSubset(
    region=instance .nodeSets['INTERFACE-NODES SIDE1'],
    position=_visualization .ELEMENT_NODAL
).values

T1 = _np .array([ (0, value .data[1], value .data[3]) for value in values]) [ranks]
T1 [: ,0] = crack_advance

with open('stress_vector-side1.dat', 'w') as fout:
    fout.write('#crack_advance[mm] Tn[MPa] Tt[MPa]
')
 np savetxt (fout, T1)

# Save nodal data in ceramic side
values = field.getSubset(
    region=instance .nodeSets['INTERFACE-NODES SIDE2'],
    position=_visualization .ELEMENT_NODAL
).values

T2 = _np .array([ (0, value .data[1], value .data[3]) for value in values]) [ranks]
T2 [: ,0] = crack_advance

with open('stress_vector-side2.dat', 'w') as fout:
    fout.write('#crack_advance[mm] Tn[MPa] Tt[MPa]
')
 np savetxt (fout, T2)

odb.close()
C Source code of 3D DBC simulation

1 from abaqus import mdb as _mdb
2 import part as _part
3 import mesh as _mesh
4 import section as _section
5 import material as _material
6 import step as _step
7
8
9
10 class ceramic:
11     length = 0.750 # mm
12     thickness = 0.317 # mm
13     young = 330000. #GPa
14     poisson = 0.22 # -
15     cte = 8.0e-6 # /dC
16
17 class copper:
18     length = 0.650 # mm
19     thickness = 0.127 # mm
20     young = 127000. # MPa
21     poisson = 0.33 # -
22     cte = 16.5e-6 # /dC
23     R0 = 43.2 # MPa
24     C1 = 27000. # MPa
25     D1 = 1850. # -
26     C2 = 1120. # MPa
27     D2 = 14.8 # -
28     Q = 30. # MPa
29     b = 12.8 # -
30
31 class interface:
32     k = 1.e8 # N/mm
33     Tmax = 150. # MPa
34     Gc = 0.115 # N/mm
35     viscosity = 0.001
36
37 class mesh:
38     size = 0.010 # mm
39
40 class loads:
41     class cooling:
42         initial = 1000.0 # dC
43         final = 0.0 # dC
44         duration = 1.0 # s
45
46     class cycles:
47         mininal = -50.0 # dC
48         maximal = 250.0 # dC
49         duration = 1.0 # s (per cycle)
50         count = 100 # cycles
51
52 name = 'simulation'
53 _model_ = None
54
55
56
```python
def create():
    ''''
    Create the simulation model in memory.
''''

This function do not take any parameter but uses the module parameters
classes to generate the model.

Example
-------

>>> import model
>>> model.mesh.size = 0.01 # Change some parameters
>>> model.create()

Warning
-------

The model is regenerated at each call using current module parameters
''''

# Delete previous generation and create a new one

global _model_
if _model_ is not None:
    _mdb.models.changeKey(fromName=_model_.name, toName='To-Delete')
model = _mdb.Model(name=name)

del _mdb.models['Model-1' if _model_ is None else 'To-Delete']

# Create the geometry
# - generate the ceramic part
sketch = model.ConstrainedSketch(name='sketch', sheetSize=1.0)
sketchrectangle(
    point1=(0., 0.),
    point2=(ceramic.length, ceramic.thickness)
)

sample = model.Part(
    name='Sample',
    dimensionality=_part.threeD,
    type=_part.DEFORMABLE_BODY
)
sample.BaseSolidExtrude(sketch=sketch, depth=ceramic.length)
del model.sketches['sketch']

# - generate the copper part
transform = sample.MakeSketchTransform(
    sketchPlane=sample faces[1],
    sketchUpEdge=sample edges[4],
    sketchPlaneSide=_part.SIDE1,
    sketchOrientation=_part RIGHT,
    origin=(0., ceramic.thickness, 0.)
)
sketch = model.ConstrainedSketch(
    name='sketch',
    sheetSize=1.0,
    transform=transform
)
sketchrectangle(
    point1=(0., 0.),
    point2=(copper.length, copper.length)
)
```
```csharp
116 sample.SolidExtrude(
117     sketchPlane=sample faces[1],
118     sketchUpEdge=sample edges[4],
119     sketchPlaneSide=_part.SIDE1,
120     sketchOrientation=_part RIGHT,
121     sketch=sketch,
122     depth=copper.thickness
123 )
124
125 # - create partitions
126 sample.PartitionCellByExtrudeEdge(
127     line=sample edges[2],
128     cells=sample cells,
129     edges=(sample edges[6], ),
130     sense=_part.REVERSE
131 )
132 sample.PartitionCellByExtendFace(
133     extendFace=sample faces[4],
134     cells=sample cells.getSequenceFromMask mask='[#1]', )
135 )
136 sample.PartitionCellByExtendFace(
137     extendFace=sample faces[8],
138     cells=sample cells.getSequenceFromMask mask='[#2]', )
139 )
140 sample.PartitionCellByExtendFace(
141     extendFace=sample faces[12],
142     cells=sample cells.getSequenceFromMask mask='[#2]', )
143 )
144
145
146 # Create the assembly
147 assembly = model rootAssembly
148 instance = assembly Instance(name='Sample', part=sample, dependent=True)
149
150
151 # Mesh the volume
152 elem_type = _mesh.ElemType(elemCode=_mesh.C3D8, elemLibrary=_mesh.STANDARD)
153 sample.set(sample cells,) elemTypes=(elem_type,))
154 sample.setMeshControls(regions=sample cells, elemShape=_mesh.QUAD)
155
156 sample.seedPart(size=mesh.size)
157 sample generateMesh()
158
159 # Create the interface mesh
160 elem_count = len(sample_elements)
161 elem_type = _mesh.ElemType(
162     elemCode=_mesh.COH3D8,
163     elemLibrary=_mesh.STANDARD,
164     viscosity=interface-viscosity
165 )
166 region = _parttalion(
167     side1Faces=sample faces.getSequenceFromMask mask='[#2000]', )
168 )
169 sample.insertElements(faces=region)
170
171 interface_set = sample. Set(
172     name='INTERFACE',
173     elements=sample elements[elem_count:]
174 )
```
175 sample.set( regions=interface_set, elemTypes=(elem_type,) )

176

177

178 # Create sets

179 bulk_set = sample.Set(name='BULK', nodes=sample.nodes)

180 ceramic_set = sample.Set(name='CERAMIC', cells=sample cells[:4])

181 copper_set = sample.Set(name='COPPER', cells=sample cells[4:])

182

183 get_faces = sample.faces.getSequenceFromMask

184 sample.Set(name='X0', faces =get_faces (mask=([ #84800 ], ), ))

185 sample.Set(name='Y0', faces =get_faces (mask=([ #400224 ], ), ))

186 sample.Set(name='Z0', faces =get_faces (mask=([ #1008080 ], ), ))

187

188

189 # Assign materials

190 # - ceramic elastic behavior

191 mat =model.Material(name='CERAMIC')

192 mat.Elastic( table=((ceramic.young, ceramic.poisson), ) )

193 mat.Expansion( table=((ceramic.cte, ), ))

194 model.HomogeneousSolidSection( name='Ceramic-Section', material='CERAMIC')

195 sample SectionAssignment( region=ceramic_set, sectionName='Ceramic-Section')

196

197 # - copper elasto-plastic behavior

198 mat =model.Material(name='COPPER')

199 mat.Elastic( table=((copper.young, copper.poisson), ) )

200 mat.Expansion( table=((copper.cte, ), ))

201 mat.Plastic(

202 hardening=_material.COMBINED,
203 dataType=_material.PARAMETERS,
204 numBackstresses=2,
205 table=((copper.R0, copper.C1, copper.D1, copper.C2, copper.D2), )
206 )

207 mat.plastic.CyclicHardening(
208 parameters=True,
209 table=((copper.R0, copper.Q, copper.b), )
210 )

211 model.HomogeneousSolidSection( name='Copper-Section', material='COPPER')

212 sample SectionAssignment( region=copper_set, sectionName='Copper-Section')

213

214 # - interface behavior

215 czm =model.Material(name='INTERFACE')

216 czm.Elastic(
217 table=((interface.k, interface.k, interface.k), ),
218 type=_material.TRACTION
219 )

220 czm.QuadsDamageInitiation(
221 table=((interface.Tmax, interface.Tmax, interface.Tmax), )
222 )

223 czm.quadsDamageInitiation.DamageEvolution(
224 table=((interface.Gc, ), ),
225 type=_material.ENERGY
226 )

227 model.CohesiveSection(
228 name='Int-Section',
229 material='INTERFACE',
230 response=_section.TRACTION_SEPARATOR
231 )

232 sample SectionAssignment( region=interface_set, sectionName='Int-Section')

233
234

# Create simulation steps

235

236

# - first cooling

237

step = modelSlowSlowStep(

238

name = 'COOLING',

239

nlgeom = True,

240

previous='Initial',

241

initialInc=loads.cooling duration/20.,

242

maxInc=loads.cooling duration/20.,

243

maxNumInc=int(1e6),

244

minInc=loads.cooling duration/1e30,

245

timePeriod=loads.cooling duration

246

)

247

step control.setValues(

248

allowPropagation = True,

249

resetDefaultValues = False,

250

timeIncrementation =(4.0, 8.0, 9.0, 16.0, 10.0, 4.0, 12.0, 50.0, 6.0, 3.0, 50.0)

251

)

252

model .fieldOutputRequests [ "F-Output-1 "] .setValues(

253

variables=( 'U', 'E', 'PE', 'PEEQ', 'S', 'SDEG', 'STATUS'),

254

timeInterval=loads.cooling duration/10.

255

)

256

# - cycles

257

258

kwargs = dict(

259

nlgeom = True,

260

initialInc=loads .cycles .duration/10.,

261

maxInc=loads .cycles .duration/10.,

262

maxNumInc = int(1e6),

263

minInc=loads .cycles .duration/1e30,

264

timePeriod=loads .cycles .duration/2.

265

)

266

previous = 'COOLING'

267

268

for i in range(1, loads .cycles .count+1):

269

step ='CYCLE'+str(i)

270

model .StaticStep( name=step+' _MIN', previous=previous, **kwargs)

271

model .StaticStep( name=step+' _MAX', previous=step+' _MIN', **kwargs)

272

previous = step+' _MAX'

273

274

if loads .cycles .count:

275

276

model .fieldOutputRequests ['F-Output-1'] .setValuesInStep(

277

stepName='CYCLE1_MIN',

278

timeInterval=loads .cycles .duration/2.

279

)

280

model .historyOutputRequests ['H-Output-1'] .setValuesInStep(

281

stepName='CYCLE1_MIN',

282

timeInterval=loads .cycles .duration/2.

283

)

284

285

# Create Boundary Conditions

286

287

model.DisplacementBC(

288

name='U1',

289

createStepName ='Initial',

290

region=instance .sets ['X0'],

291

u1=0.0

292

)
```python
293     model.DisplacementBC(
294         name='U2',
295         createStepName='Initial',
296         region=instance.sets['YO'],
297         u2=0.0
298     )
299     model.DisplacementBC(
300         name='U3',
301         createStepName='Initial',
302         region=instance.sets['ZO'],
303         u3=0.0
304     )
305
306 # Apply temperature
307 field = model.Temperature(
308     name='TEMP',
309     createStepName='Initial',
310     region=instance.sets['BULK'],
311     magnitudes=(loads.cooling(initial, )
312 )
313 field.setValuesInStep(
314     stepName='COOLING',
315     magnitudes=(loads.cooling(final, )
316 )
317 for i in range(1, loads.cycles.count+1):
318     field.setValuesInStep(
319         stepName='CYCLE'+str(i)+'_MIN',
320         magnitudes=(loads.cycles.minimal, )
321     )
322     field.setValuesInStep(
323         stepName='CYCLE'+str(i)+'_MAX',
324         magnitudes=(loads.cycles.maximal, )
325     )
326
327 # Store the model globaly
328 __model__ = model
329
330
def write_input():
331     '''
332     Write the simulation inp file from the created model
333     '''
334     Example
335     -------
336     >>> import model
337     >>> model.create()
338     >>> model.write_input()
339
340
341 Raises
342     -------
343     RuntimeError if the model has not been generated
344     '''
345     if __model__ is None:
346         raise RuntimeError('Model not generated')
347
348     job = _mdb.Job(name=__model__.name, model=__model__)
349     job.writeInput(consistencyChecking=False)
350 del _mdb.jobs[__model__.name]
```
References

[Ben Kaabar, 2015] Ben Kaabar, A. (2015). *Durabilité des assemblages céramique-métal employés en électronique de puissance*. Thesis, INSA de Lyon.

[Dugdale, 1960] Dugdale, D. S. (1960). Yielding of steel sheets containing slits. *Journal of The Mechanics and Physics of Solids*, 8:100–104.

[Gaiser et al., 2015] Gaiser, P., Klingler, M., and Wilde, J. (2015). Fracture mechanical modeling for the stress analysis ofdbc ceramics. *2015 16th International Conference on Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Microsystems*, pages 1–6.