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.
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")
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