Example 6 - 3D model
Problem description
Figure 1 depicts the model considered in this example. We are interested in the gradient of the strain energy with respect to the position of the cylindrical hole. Configurational forces are used to estimate this energy gradient. The configurational forces are computed using ConForce. To validate this approach, the energy gradient is estimated numerically by a difference quotient, for which two additional FEM simulations with slightly displaced holes are run.
Figure 1: 3D model with a cylindrical hole. The left side is fixed. A displacement \(u_{y}\) is applied on the right side. The regions with radii \(r_{1}\) and \(r_{2}\) are used to compute the configurational forces. The dimensions are listed in section Geometric dimensions.
The mesh is shown in Figure 2. We use 20-node bricks with reduced integration. The image shows a red-, green-, and cream-colored region. The cream-colored region has the radius \(r_{1}\) and the green region has the radius \(r_{2}\). The configurational forces are computed for all nodes inside two domains:
The smaller domain \(\mathcal{A}\) contains all nodes of the cream-colored region.
The larger domain \(\mathcal{B}\) contains all nodes of the green- and cream-colored region.
Figure 2: The FEM mesh consists of 20-node brick elements with reduced integration (C3D20R).
Simulation
>>> import os
>>> import json
>>> import subprocess
>>> import numpy as np
>>> import scipy.optimize as optimize
>>> HOME_DIR = os.path.abspath(".")
>>> os.chdir(__file__ + "/..")
The model is automatically build and simulated by an Abaqus script. We provide the name of the model and the position of the hole by cx and cy.
>>> model_name = "bending_model_0"
>>> cx = 50 # mm
>>> cy = 40 # mm
>>> subprocess.call(f"abaqus cae noGui=example_6_abaqus_script.py -- {model_name} {cx:.3f} {cy:.3f}", shell=True)
0
The Abaqus script writes results to a json file. The command below loads the results.
>>> with open(f"{model_name}_results.json", "r", encoding="UTF-8") as fh:
... results = json.load(fh)
Material
The compressible Neo-Hookean material model for hyperelastic materials is used. The first parameter
>>> results["C10"] # MPa
86.0
defines the deviatoric behavior. The second parameter
>>> results["D1"] # 1/MPa
0.0012
controls the behavior under a hydrostatic load.
Geometric dimensions
The dimensions of the model shown in Figure 1 are:
>>> results["lx"] # mm
100.0
>>> results["ly"] # mm
50.0
>>> results["lz"] # mm
5.0
>>> results["cx"] # mm
50.0
>>> results["cy"] # mm
40.0
>>> results["r"] # mm
5.0
>>> results["r1"] # mm
6.0
>>> results["r2"] # mm
8.0
Applied displacement
>>> results["uy"] # mm
10.0
Configurational forces
The configurational forces in x- and y-direction are computed for two domains. Domain \(\mathcal{A}\) has configurational forces of
>>> CF_A = results["CF_A"]
>>> # in mJ/mm**2
>>> CF_A[:-1]
[4.903..., -16.89...]
Domain \(\mathcal{B}\) is larger and configurational forces of
>>> CF_B = results["CF_B"]
>>> # in mJ/mm**2
>>> CF_B[:-1]
[4.892..., -16.90...]
matches the results of \(\mathcal{A}\) with less than 0.3% deviation. We consider only the volume integral (23) and neglect the surface influence on the configurational forces.
Numerical energy gradient
The energy gradient with respect to the position of the hole can also be found by a numerical derivation. For the numerical derivation, we move the hole once by a small dx
>>> model_name_dx = "bending_model_dx"
>>> dx = 0.1 # mm
>>> cx_dx = cx + dx
>>> subprocess.call(f"abaqus cae noGui=example_6_abaqus_script.py -- {model_name_dx} {cx_dx:.3f} {cy:.3f}", shell=True)
0
and once by a small dy:
>>> model_name_dy = "bending_model_dy"
>>> dy = 0.1 # mm
>>> cy_dy = cy + dy
>>> subprocess.call(f"abaqus cae noGui=example_6_abaqus_script.py -- {model_name_dy} {cx:.3f} {cy_dy:.3f}", shell=True)
0
The simulation computes the strain energy (ALLSE) and writes it into json files.
>>> with open(f"{model_name_dx}_results.json", "r", encoding="UTF-8") as fh:
... results_dx = json.load(fh)
>>> with open(f"{model_name_dy}_results.json", "r", encoding="UTF-8") as fh:
... results_dy = json.load(fh)
Next, the energies (ALLSE) are read from the json files.
>>> ALLSE_0 = results["ALLSE"]
>>> ALLSE_dx = results_dx["ALLSE"]
>>> ALLSE_dy = results_dy["ALLSE"]
With the energies, we can use a forward difference quotient to numerically estimate the energy gradient.
>>> Gx = (ALLSE_dx - ALLSE_0) / dx
>>> Gy = (ALLSE_dy - ALLSE_0) / dy
>>> # in mJ/mm**2
>>> [Gx, Gy]
[4.887..., -17.1...]
Conclusion
The energy gradient with respect to x computed by the configurational forces is almost identical to the result of the numerical derivation. The energy gradient with respect to y is also in a good agreement, although the configurational forces method is slightly smaller than the result obtained by the numerical derivation.
To minimize the strain energy, the hole should be moved in the opposite direction of the gradient. Consequently, moving the hole up and to the left would decrease the strain energy. This seems plausible because the bending stiffness (and the energy) will be less if the hole is at the top compared to the hole being in the middle. Furthermore, the bending moment (and the energy density) is greater on the left side. A hole on the left side reduces the total energy more than a hole on the right side.
Abaqus script
from __future__ import print_function
import sys
import os
import json
import numpy as np
import abaqus as abq
import abaqusConstants as abqConst
import caeModules as cae
# 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
# parameters
C10 = 86.
D1 = 1.2e-3
def main(
name="bending_beam_model",
lx=100., ly=50., lz=5.,
cx=50., cy=40., r=5.,
r1=6., r2=8.,
uy=10.
):
"""
creates a model with a cylindrical hole that is fixed on the left side.
A displacement is applied on the right side.
:param name: name of the model, the job and odb files
:param lx: length of the model
:param ly: height of the model
:param lz: thickness of the model
:param cx: x-position of the cylindrical hole
:param cy: y-position of the cylindrical hole
:param r: radius of the cylindrical hole
:param r1: radius (r1 > r) of region "1" around the cylindrical hole for the evaluation of computational forces
:param r2: radius (r2 > r1) of region "2" around the cylindrical hole for the evaluation of computational forces
:param uy: applied displacement on the right side of the model
:return: dictionary containing the strain energy (ALLSE) and the configurational forces of region "1" and "2" (CF_1 and CF_2)
"""
# put all results in this dictionary
results = dict(
name=name,
lx=lx,
ly=ly,
lz=lz,
cx=cx,
cy=cy,
r=r,
r1=r1,
r2=r2,
uy=uy,
C10=C10,
D1=D1
)
# Clear model database
abq.Mdb()
# Create model
model = abq.mdb.Model(name=name)
del abq.mdb.models["Model-1"]
# create part
part = model.Part(
name='Part-1',
dimensionality=abqConst.THREE_D,
type=abqConst.DEFORMABLE_BODY
)
sketch = model.ConstrainedSketch(name='sketch', sheetSize=lx)
sketch.rectangle(point1=(0.0, 0.0), point2=(lx, ly))
sketch.CircleByCenterPerimeter(center=(cx, cy), point1=(cx+r, cy))
part.BaseSolidExtrude(sketch=sketch, depth=lz)
# partition part
partition_sketch = model.ConstrainedSketch(name='partition', sheetSize=lx)
partition_sketch.CircleByCenterPerimeter(center=(cx, cy), point1=(cx+r1, cy))
partition_sketch.CircleByCenterPerimeter(center=(cx, cy), point1=(cx+r2, cy))
part.PartitionFaceBySketch(
faces=part.faces.findAt(((cx, cy*0.2, 0), )),
sketch=partition_sketch
)
part.PartitionCellByExtrudeEdge(
cells=part.cells[:],
edges=part.edges.findAt(
((cx, cy+r1, 0), ),
((cx, cy+r2, 0), ),
),
line=part.edges.findAt((0., 0., lz*0.2), ),
sense=abqConst.FORWARD
)
# Set
part_all_set = part.Set(name="all", cells=part.cells)
# create material, section, section assignments
material = model.Material(name="material")
material.Hyperelastic(
type=abqConst.NEO_HOOKE,
materialType=abqConst.ISOTROPIC,
testData=abqConst.OFF,
volumetricResponse=abqConst.VOLUMETRIC_DATA,
table=((C10, D1), ))
section = model.HomogeneousSolidSection(
name="section",
material="material",
thickness=None
)
part.SectionAssignment(
region=part_all_set,
sectionName=section.name,
)
# mesh part
part.seedPart(size=lx/50, deviationFactor=0.1, minSizeFactor=0.1)
num_nodes_circle = int((2*r2*np.pi) / (lx/100))
part.seedEdgeByNumber(
edges=part.edges.findAt(((cx, cy+r, 0), ),),
number=num_nodes_circle
)
part.seedEdgeByNumber(
edges=part.edges.findAt(((cx, cy+r1, 0), ),),
number=num_nodes_circle
)
part.seedEdgeByNumber(
edges=part.edges.findAt(((cx, cy+r2, 0), ),),
number=num_nodes_circle
)
part.setMeshControls(
regions=part.cells.findAt(((cx, cy + (r + r1) / 2, lz*0.2),)),
algorithm=abqConst.MEDIAL_AXIS
)
part.setElementType(
regions=part_all_set,
elemTypes=(
cae.mesh.ElemType(elemCode=abqConst.C3D20R, elemLibrary=abqConst.STANDARD),
cae.mesh.ElemType(elemCode=abqConst.C3D15, elemLibrary=abqConst.STANDARD),
cae.mesh.ElemType(elemCode=abqConst.C3D10, elemLibrary=abqConst.STANDARD)
)
)
part.generateMesh()
# create instance
assembly = model.rootAssembly
instance = assembly.Instance(name='Part-1-1', part=part, dependent=abqConst.ON)
rp_feature = assembly.ReferencePoint(point=(lx, ly/2, lz/2))
rp = assembly.referencePoints[rp_feature.id]
# create sets
assembly.Set(
name="REGION_A",
cells=instance.cells.findAt(((cx, cy + (r + r1) / 2, lz*0.2),))
)
assembly.Set(
name="REGION_B",
cells=instance.cells.findAt(
((cx, cy + (r + r1) / 2, lz*0.2),),
((cx, cy + (r1 + r2) / 2, lz*0.2),),
)
)
node_rp_set = assembly.Set(
name="node_rp",
referencePoints=(rp, )
)
node_000_set = assembly.Set(
name="node_000",
vertices=instance.vertices.findAt(((0., 0., 0.),))
)
face_0_set = assembly.Set(
name="face_0",
faces=instance.faces.findAt(((0., ly*0.2, lz*0.2),))
)
face_1_set = assembly.Set(
name="face_1",
faces=instance.faces.findAt(((lx, ly*0.2, lz*0.2),))
)
surface_1 = assembly.Surface(
name='face_1',
side1Faces=face_1_set.faces
)
# coupling
model.Coupling(
name='Constraint-1',
surface=surface_1,
controlPoint=node_rp_set,
influenceRadius=abqConst.WHOLE_SURFACE,
couplingType=abqConst.KINEMATIC,
u1=abqConst.ON, u2=abqConst.ON, u3=abqConst.ON,
ur1=abqConst.ON, ur2=abqConst.ON, ur3=abqConst.ON
)
# create step
step = model.StaticStep(name='Step-1', previous='Initial', nlgeom=abqConst.ON)
model.fieldOutputRequests['F-Output-1'].setValues(variables=('S', 'LE', 'U', 'ENER'))
# create boundary conditions
model.PinnedBC(name='pinned', createStepName='Initial', region=face_0_set, localCsys=None)
model.DisplacementBC(name='BC', createStepName=step.name, region=node_rp_set, u2=-uy)
# save cae
abq.mdb.saveAs(pathName=name + ".cae")
#################################
# SIMULATION AND POSTPROCESSING #
#################################
# create a job, start the simulation and wait until the simulation completed
job = abq.mdb.Job(
name=name,
model=model.name,
nodalOutputPrecision=abqConst.FULL
)
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",
method="mbf",
e_expression="SENER"
)
# compute results
# configurational forces
fo_CF = odb.steps['Step-1'].frames[-1].fieldOutputs["CONF_FORCE"]
results["CF_A"] = np.sum([
value.data
for value in fo_CF.getSubset(region=odb.rootAssembly.nodeSets["REGION_A"]).values
], axis=0).tolist()
results["CF_B"] = np.sum([
value.data
for value in fo_CF.getSubset(region=odb.rootAssembly.nodeSets["REGION_B"]).values
], axis=0).tolist()
# strain energy
results["ALLSE"] = odb.steps['Step-1'].historyRegions['Assembly ASSEMBLY'].historyOutputs['ALLSE'].data[-1][-1]
return results
if __name__ == '__main__':
name = str(sys.argv[-3])
cx = float(sys.argv[-2])
cy = float(sys.argv[-1])
result = main(
name=name,
cx=cx,
cy=cy
)
with open(name + "_results.json", "w") as fh:
json.dump(result, fh, indent=1)