Example 2 - CT-specimen with linear elastic material behaviour

Problem description

Kolednik [1] suggests the following model setup for a compact tension (CT) test. Contrary to Kolednik, in this example a linear elastic material behaviour is used instead of an elastic plastic behaviour. This allows a comparison with Anderson’s theoretical prediction for CT specimens [2]. The regions A and B have a structured mesh with a mesh size of .15 mm, and 0.167 mm. Region C and D have an approximate mesh size of 1 mm. Region D is defined for Example 3 where elastic-plastic material behaviour is investigated.

scheme

Model of a compact tension specimen

Geometric dimensions:

>>> w = 50  # mm
>>> a = 27  # mm
>>> b = 25  # mm
>>> h = 60  # mm

Material properties:

>>> nu = 0.3
>>> E = 200_000 # MPa

Applied displacement:

>>> u = 0.5  # mm

Simulation

A *.inp file of the model is provided in the example folder. The model evaluation is scripted and there is no need to open the ConForce-plugin manually. First, change to the directory where the *.inp file is located.

>>> import os
>>> import numpy as np
>>> HOME_DIR = os.path.abspath(".")
>>> os.chdir(__file__ + "/..")

Next, call the Abaqus script (example_2_abaqus_script). The script simulates the *.inp file and writes a results.json file.

>>> import subprocess
>>> subprocess.call("abaqus cae noGUI=example_2_abaqus_script.py", shell=True)
0

After the simulation, the result file is loaded to fetch the reaction force load. This force is needed to reach the defined displacement u.

>>> import json
>>> with open("results.json", "r", encoding="UTF-8") as fh:
...     results = json.load(fh)
>>> load = results["reaction_force"]
>>> np.around(load, 3)  # N
56319.531

Theory

Anderson [2] defines an analytical solution for this problem.

>>> geometry_factor = (
...     (2 + a/w)
...     * (
...         0.866
...         + 4.64 * (a/w)
...         - 13.32 * (a/w)**2
...         + 14.72 * (a/w)**3
...         - 5.6 * (a/w)**4
...     )
...     / (1 - a/w)**(3./2)
... )
>>> geometry_factor
10.821389667870234

The stres intensity factor k is:

>>> k = (
...     load  # N
...     / (b * w**0.5)  # 1/mm**(3/2)
...      * geometry_factor
... )  # N / mm**(3/2) = MPa mm**0.5
>>> np.around(k, 3)  # MPa mm**0.5
3447.601

For a plane strain state, the J-Integral is computed from k as:

>>> J_theory = k**2 * (1 - nu**2) / E
>>> np.around(J_theory, 3)  # mJ/mm**2
54.081

Abaqus J-Integral

Abaqus computes the J-Integral using the virtual crack extension method by Parks [3]. Abaqus computes the J-Integral over several contours. The regions are defined in the image and correspond to the following contour indices

  • region A corresponds to the contour index 21

  • region B corresponds to the contour index 57

According to Abaqus the J-Integral for region A is

>>> J_in_a_abaqus = results["J"]["J at J_NEAR_CRACK_TIP_Contour_21"]
>>> np.around(J_in_a_abaqus, 3)  # mJ/mm**2
54.503

for region B

>>> J_in_b_abaqus = results["J"]["J at J_NEAR_CRACK_TIP_Contour_57"]
>>> np.around(J_in_b_abaqus, 3)  # mJ/mm**2
54.635

and for region FAR_FIELD

>>> J_in_far_abaqus = results["J"]["J at J_FAR_FAR_FIELD_Contour_1"]
>>> np.around(J_in_far_abaqus, 3)  # mJ/mm**2
54.765

This is in good aggreement with the prediction made by Anderson. Furthermore, the path-independency of the J-integral is fulfilled well.

Configurational forces

ConForce can predict J-Integrals by summing up configurational forces inside a region. The resulting configurational force is projected onto the crack extension direction. In this case the crack extension direction is simply the x-axis.

The configurational forces for the regions A, B, and FAR_FIELD show good aggreement with Anderson and Abaqus.

Note

The configurational forces have a negative sign.

>>> (cfx_in_a, _) = results["CF"]["A"]
>>> np.around(cfx_in_a, 3)  # mJ/mm**2
-54.485
>>> (cfx_in_b, _) = results["CF"]["B"]
>>> np.around(cfx_in_b, 3)  # mJ/mm**2
-54.633
>>> (cfx_in_far, _) = results["CF"]["FAR_FIELD"]
>>> np.around(cfx_in_far, 3)  # mJ/mm**2
-54.808

Like conventional forces, the configurational forces fullfill the equilibrium of forces. The configuraional forces summed up for the entire model are zero.

>>> (cfx_whole_model, _) = results["CF"]["ALL_NODES"]
>>> np.around(cfx_whole_model, 3)
0.0

Comparison

The figure compares the J-Integral of Abaqus and Anderson with the negative configurational forces of ConForce. Abaqus and ConForce show a good aggreement for all contours in the regions A and B.

>>> fig = compare_J_and_negative_cfx(results, J_theory)
>>> save_fig(fig, "example_2_images/01_comparison_over_contours.png")
comparison

References

Change to home directory

>>> os.chdir(HOME_DIR)

Abaqus script

"""
This is an Abaqus script for the example in :py:mod:`example_2_CT_specimen_linear_elastic`.

To simulated and write a result.json file open a shell and type:

.. code-block:: console

    abaqus cae noGUI="example_2_abaqus_script.py"

"""
from __future__ import print_function

import sys
import os
import json

import numpy as np

import abaqus as abq

# append path to the folder where conforce_abq is located and import conforce_abq afterward
sys.path.append(os.path.abspath("../.."))
from conforce_abq.main import apply


def simulate_and_save_results(inp_file_path="CT_specimen_CPE4.inp"):
    """
    This function:
    - simulate the input file
    - computes the configurational forces in the odb
    - writes results to a file called "results.json"

    :param inp_file_path: str, path to input file to simulate
    :return: updated Odb with configurational forces as field output
    """

    # create a job from the input file, start the simulation and wait until the simulation completed
    job_name = os.path.basename(inp_file_path).split(".")[0]
    job = abq.mdb.JobFromInputFile(
        name=job_name,
        inputFileName=inp_file_path
    )
    job.submit()
    job.waitForCompletion()

    # open odb with **readOnly=False**
    odb = abq.session.openOdb(job_name + ".odb", readOnly=False)

    # apply conforce plugin and request configurational forces as field output
    odb = apply(
        odb,
        request_CF=True,
        CF_name="CONF_FORCE_LIN_EL",
        method="mbf",
        e_expression="SENER"
    )

    # put all results in this dictionary
    results = dict()

    # extract CF field output
    fo_CF = odb.steps['Loading'].frames[-1].fieldOutputs["CONF_FORCE_LIN_EL"]

    results["CF"] = dict()
    for set_name, region in odb.rootAssembly.nodeSets.items():
        fo_CF_in_region = fo_CF.getSubset(
            region=region
        )

        results["CF"][set_name] = np.sum([
            value.data
            for value in fo_CF_in_region.values
        ], axis=0).tolist()

    # extract history outputs
    (ho_assembly, ho_J_integral, ho_bc_lower, ho_bc_upper) = [
        history_region.historyOutputs
        for history_region in odb.steps['Loading'].historyRegions.values()
    ]

    # extract J-Integrals at the last frame
    results["J"] = {
        key: value.data[-1][-1]
        for key, value in ho_J_integral.items()
    }

    # extract reaction forces
    results["reaction_force"] = ho_bc_upper['RF2'].data[-1][-1]

    # save results
    with open("results.json", "w") as fh:
        json.dump(results, fh, indent=4)

    return odb


if __name__ == '__main__':
    # delegate output to console
    stdout_default, stderr_default = sys.stdout, sys.stderr
    sys.stdout, sys.stderr = sys.__stdout__, sys.__stderr__
    try:
        # simulate
        simulate_and_save_results()

    finally:
        # restore default output streams
        sys.stdout, sys.stderr = stdout_default, stderr_default