Note

This notebook can be downloaded here: indentation_2D_Fibre.ipynb

2D indentation of a fibre in a matrix

%load_ext autoreload
%autoreload 2
import matplotlib.pyplot as plt
import matplotlib as mpl
import hardness as hd
import argiope as ag
import pandas as pd
import numpy as np
import os, subprocess, time, local_settings, time
%matplotlib nbagg

mpl.rcParams['grid.color'] = 'k'
mpl.rcParams['grid.linestyle'] = ':'
mpl.rcParams['grid.linewidth'] = 0.5
mpl.rcParams['contour.negative_linestyle'] = 'solid'

# USEFUL FUNCTIONS
def create_dir(path):
  try:
    os.mkdir(path)
  except:
    pass

Settings

# SETTINGS
workdir   = "workdir/"
outputdir = "outputs/"
label   = "indentation_2D"

create_dir(workdir)
create_dir(workdir + outputdir)

Model definition

#-------------------------------------------------------------------------------
# MESH DEFINITIONS
def element_map(mesh):
    mesh.elements.loc[mesh.elements.type.argiope == "tri3", ("type", "solver", "")] = "CAX3"
    mesh.elements.loc[mesh.elements.type.argiope == "quad4", ("type", "solver", "")] = "CAX4"
    return mesh


def sample_material_map(mesh):
    mesh.elements.loc[mesh.elements.sets.FIBRE, "materials"] = "FIBRE_MAT"
    mesh.elements.loc[mesh.elements.sets.MATRIX, "materials"] = "MATRIX_MAT"
    return mesh

def indenter_material_map(mesh):
    mesh.elements["materials"] = "INDENTER_MAT"
    return mesh


parts = {
    "sample" : hd.models.SampleFibre2D(Rf = 1.,
                                   ly1 = 1., ly2 = 10.,
                                   Nx = 16, Ny = 8,
                                   Nr = 8, Nt = None,
                                   gmsh_path = "gmsh",
                                   file_name = "dummy",
                                   workdir = workdir,
                                   gmsh_space = 2,
                                   gmsh_options = "-algo 'delquad'",
                                   element_map = element_map,
                                   material_map = sample_material_map),

    "indenter" : hd.models.SpheroconicalIndenter2D(
                                   R = 1.,
                                   psi= 70.3,
                                   r1 = 1.,
                                   r2 = 3.,
                                   r3 = 3.,
                                   lc1 = .1,
                                   lc2 = .5,
                                   rigid = False,
                                   gmsh_path = "gmsh",
                                   file_name = "dummy",
                                   workdir = workdir,
                                   gmsh_space = 2,
                                   gmsh_options = "-algo 'delquad'",
                                   element_map = element_map,
                                   material_map = indenter_material_map)}

materials = [ag.materials.ElasticPerfectlyPlastic(
                                  label = "MATRIX_MAT",
                                  young_modulus = .1,
                                  poisson_ratio = .3,
                                  yield_stress = .001),
             ag.materials.Elastic(label = "INDENTER_MAT",
                                  young_modulus = 1.,
                                  poisson_ratio = .3),
             ag.materials.Elastic(label = "FIBRE_MAT",
                                  young_modulus = .4,
                                  poisson_ratio = .3)]

#-------------------------------------------------------------------------------
# STEP DEFINTIONS
steps = [
        hd.models.Step2D(name = "LOADING1",
                         control_type = "disp",
                         duration = 1.,
                         kind = "adaptative",
                         nframes = 100,
                         controlled_value = -0.2,
                         field_output_frequency = 99999),
        hd.models.Step2D(name = "UNLOADING1",
                         control_type = "force",
                         duration = 1.,
                         kind = "adaptative",
                         nframes = 100,
                         controlled_value = 0.,
                         field_output_frequency = 99999),
        hd.models.Step2D(name = "RELOADING1",
                         control_type = "disp",
                         duration = 1.,
                         kind = "adaptative",
                         nframes = 100,
                         controlled_value = -0.2,
                         field_output_frequency = 99999),
        hd.models.Step2D(name = "LOADING2",
                         control_type = "disp",
                         duration = 1.,
                         kind = "adaptative",
                         nframes = 100,
                         controlled_value = -0.4,
                         field_output_frequency = 99999),
        hd.models.Step2D(name = "UNLOADING2",
                         control_type = "force",
                         kind = "adaptative",
                         duration = 1.,
                         nframes = 50,
                         controlled_value = 0.,
                         field_output_frequency = 99999)
        ]

model0 = hd.models.Indentation2D(label = label,
                      parts = parts,
                      steps = steps,
                      materials = materials,
                      solver = "abaqus",
                      solver_path = local_settings.ABAQUS_PATH,
                      workdir = workdir,
                      verbose = True)
print("1: Preprocessing ----------------------------------")
%time model0.write_input()
print("2: Processing -------------------------------------")
%time model0.run_simulation()
print("3: Postprocessing ---------------------------------")
%time model0.postproc()
print("4: Saving model -----------------------------------")
%time model0.save(workdir + "model.pcklz")
1: Preprocessing ----------------------------------
CPU times: user 1.48 s, sys: 20 ms, total: 1.5 s
Wall time: 2.46 s
2: Processing -------------------------------------
<Running "indentation_2D" using abaqus>
(b'Abaqus JOB indentation_2DnAbaqus 6.13-1nBegin Analysis Input File ProcessornSun Jan 14 22:02:21 2018nRun prenSun Jan 14 22:02:26 2018nEnd Analysis Input File ProcessornBegin Abaqus/Standard AnalysisnSun Jan 14 22:02:26 2018nRun standardnSun Jan 14 22:03:11 2018nEnd Abaqus/Standard AnalysisnAbaqus JOB indentation_2D COMPLETEDn', None)
<Ran indentation_2D: duration 54.94s>
CPU times: user 4 ms, sys: 8 ms, total: 12 ms
Wall time: 54.9 s
3: Postprocessing ---------------------------------
<Post-Processing"indentation_2D" using abaqus>
/opt/abaqus/scratch/Commands/abaqus viewer noGUI=indentation_2D_abqpp.py
(b'', None)
<Post-Processed indentation_2D: duration 11.13s>
CPU times: user 296 ms, sys: 16 ms, total: 312 ms
Wall time: 11.4 s
4: Saving model -----------------------------------
CPU times: user 88 ms, sys: 0 ns, total: 88 ms
Wall time: 89.5 ms
model = ag.utils.load(workdir + "model.pcklz")

Model checking

Mesh building and quality checking.

model.parts["indenter"].mesh.elements.head()
conn materials sets type surfaces
n0 n1 n2 n3 ALL_ELEMENTS argiope solver SURFACE
f1 f2 f3 f4
element
47 85 62 73 64 INDENTER_MAT True quad4 CAX4 False False False False
48 132 103 92 66 INDENTER_MAT True quad4 CAX4 False False False False
49 57 60 72 80 INDENTER_MAT True quad4 CAX4 False False False False
50 76 146 124 111 INDENTER_MAT True quad4 CAX4 False False False False
51 156 160 166 153 INDENTER_MAT True quad4 CAX4 False False False False
model.parts["sample"].mesh.elements.head()
conn materials sets type sets surfaces
n0 n1 n2 n3 FIBRE MATRIX argiope solver ALL_ELEMENTS SURFACE
f1 f2 f3 f4
element
68 343 359 345 0 MATRIX_MAT False True tri3 CAX3 True False False False False
69 341 344 350 0 MATRIX_MAT False True tri3 CAX3 True False False False False
70 358 365 368 0 MATRIX_MAT False True tri3 CAX3 True False False False False
71 337 352 347 0 MATRIX_MAT False True tri3 CAX3 True False False False False
72 333 351 345 0 MATRIX_MAT False True tri3 CAX3 True False False False False
parts = model.parts
i = 1
fig = plt.figure()
parts_names = parts.keys()
for name, part in parts.items():
    mesh = part.mesh
    patches = mesh.to_polycollection(edgecolor = "black", linewidth = .5, alpha = 1.)
    stats = mesh.stats()
    patches.set_array( stats.stats.max_abs_angular_deviation )
    patches.set_cmap(mpl.cm.YlOrRd)
    ax = fig.add_subplot(1, 2, i)
    ax.set_aspect("equal")
    ax.set_xlim(mesh.nodes.coords.x.min(), mesh.nodes.coords.x.max())
    ax.set_ylim(mesh.nodes.coords.y.min(), mesh.nodes.coords.y.max())
    ax.add_collection(patches)
    cbar = plt.colorbar(patches, orientation = "horizontal")
    cbar.set_label("Max Abs. Angular Deviation [$^o$]")
    plt.xlabel("$x$")
    plt.ylabel("$y$")
    plt.grid()
    plt.title(name.title())
    i+= 1
plt.show()
<IPython.core.display.Javascript object>

Simulation

Post-Processing

Time data

hist = model.data["history"]
hist.head()
Wes Wf RF Wtot dtot dtip Wei CF Wps t step F
0 0.000000e+00 0.0 -0.000000 0.000000e+00 0.000 0.000000 0.000000e+00 0.0 0.0 0.00 0 0.000000
1 1.551044e-08 0.0 -0.000025 2.913694e-08 -0.002 -0.001598 4.527072e-09 0.0 0.0 0.01 0 -0.000025
2 7.725865e-08 0.0 -0.000071 1.374515e-07 -0.004 -0.003053 2.598622e-08 0.0 0.0 0.02 0 -0.000071
3 2.503104e-07 0.0 -0.000170 3.783753e-07 -0.006 -0.004348 9.346720e-08 0.0 0.0 0.03 0 -0.000170
4 5.718712e-07 0.0 -0.000268 8.160462e-07 -0.008 -0.005641 2.087147e-07 0.0 0.0 0.04 0 -0.000268
plt.figure()
for step, group in hist.groupby("step"):
  plt.plot(-group.dtot, -group.F, label = "Step {0}".format(step))
plt.grid()
plt.legend(loc = "best")
plt.ylabel("Total force $F$, []")
plt.xlabel("Displacement, $\delta$ []")
plt.show()
<IPython.core.display.Javascript object>

Fields

model.parts["sample"].mesh.fields_metadata()
frame frame_value label part position step_label step_num
0 0 0 S I_SAMPLE node LOADING1 0
1 0 0 U I_SAMPLE node LOADING1 0
2 1 1 S I_SAMPLE node LOADING1 0
3 1 1 U I_SAMPLE node LOADING1 0
4 0 0 S I_SAMPLE node UNLOADING1 1
5 0 0 U I_SAMPLE node UNLOADING1 1
6 1 1 S I_SAMPLE node UNLOADING1 1
7 1 1 U I_SAMPLE node UNLOADING1 1
8 0 0 S I_SAMPLE node RELOADING1 2
9 0 0 U I_SAMPLE node RELOADING1 2
10 1 1 S I_SAMPLE node RELOADING1 2
11 1 1 U I_SAMPLE node RELOADING1 2
12 0 0 S I_SAMPLE node LOADING2 3
13 0 0 U I_SAMPLE node LOADING2 3
14 1 1 S I_SAMPLE node LOADING2 3
15 1 1 U I_SAMPLE node LOADING2 3
16 0 0 S I_SAMPLE node UNLOADING2 4
17 0 0 U I_SAMPLE node UNLOADING2 4
18 1 1 S I_SAMPLE node UNLOADING2 4
19 1 1 U I_SAMPLE node UNLOADING2 4
model.parts["sample"].mesh.fields_metadata()
frame frame_value label part position step_label step_num
0 0 0 S I_SAMPLE node LOADING1 0
1 0 0 U I_SAMPLE node LOADING1 0
2 1 1 S I_SAMPLE node LOADING1 0
3 1 1 U I_SAMPLE node LOADING1 0
4 0 0 S I_SAMPLE node UNLOADING1 1
5 0 0 U I_SAMPLE node UNLOADING1 1
6 1 1 S I_SAMPLE node UNLOADING1 1
7 1 1 U I_SAMPLE node UNLOADING1 1
8 0 0 S I_SAMPLE node RELOADING1 2
9 0 0 U I_SAMPLE node RELOADING1 2
10 1 1 S I_SAMPLE node RELOADING1 2
11 1 1 U I_SAMPLE node RELOADING1 2
12 0 0 S I_SAMPLE node LOADING2 3
13 0 0 U I_SAMPLE node LOADING2 3
14 1 1 S I_SAMPLE node LOADING2 3
15 1 1 U I_SAMPLE node LOADING2 3
16 0 0 S I_SAMPLE node UNLOADING2 4
17 0 0 U I_SAMPLE node UNLOADING2 4
18 1 1 S I_SAMPLE node UNLOADING2 4
19 1 1 U I_SAMPLE node UNLOADING2 4
parts = {k:part.mesh.copy() for k, part in model.parts.items() }

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_aspect("equal")
ax.set_xlim(0., 3.)
ax.set_ylim(-2., 2.)

field_num = 14
disp_num = 15
levels = np.linspace(-1.e-1, 1.e-1, 11)

for k, mesh in parts.items():
    field = mesh.fields[field_num].data.v22
    disp = mesh.fields[disp_num].data
    mesh.nodes[("coords", "x")] += disp.v1
    mesh.nodes[("coords", "y")] += disp.v2
    tri = mesh.to_triangulation()
    patches = mesh.to_polycollection(facecolor = "none",
                                     edgecolor = "black",
                                     linewidth = .5)
    grad = ax.tricontourf(tri, field, levels, cmap = mpl.cm.terrain, alpha = 1)
    ax.tricontour(tri, field, levels, colors = "white", linewidths = 1.)
    ax.add_collection(patches)
cbar = plt.colorbar(grad)
cbar.set_label("Cauchy Stress, $\sigma_{12}$")
plt.xlabel("$x$")
plt.ylabel("$y$")
#plt.grid()
<IPython.core.display.Javascript object>
<matplotlib.text.Text at 0x7fd9740942e8>
parts["indenter"].fields
[<Tensor4Field S at node; step=0 ('LOADING1'); frame=0 >,
 <Vector2Field U at node; step=0 ('LOADING1'); frame=0 >,
 <Tensor4Field S at node; step=0 ('LOADING1'); frame=1 >,
 <Vector2Field U at node; step=0 ('LOADING1'); frame=1 >,
 <Tensor4Field S at node; step=1 ('UNLOADING1'); frame=0 >,
 <Vector2Field U at node; step=1 ('UNLOADING1'); frame=0 >,
 <Tensor4Field S at node; step=1 ('UNLOADING1'); frame=1 >,
 <Vector2Field U at node; step=1 ('UNLOADING1'); frame=1 >,
 <Tensor4Field S at node; step=2 ('RELOADING1'); frame=0 >,
 <Vector2Field U at node; step=2 ('RELOADING1'); frame=0 >,
 <Tensor4Field S at node; step=2 ('RELOADING1'); frame=1 >,
 <Vector2Field U at node; step=2 ('RELOADING1'); frame=1 >,
 <Tensor4Field S at node; step=3 ('LOADING2'); frame=0 >,
 <Vector2Field U at node; step=3 ('LOADING2'); frame=0 >,
 <Tensor4Field S at node; step=3 ('LOADING2'); frame=1 >,
 <Vector2Field U at node; step=3 ('LOADING2'); frame=1 >,
 <Tensor4Field S at node; step=4 ('UNLOADING2'); frame=0 >,
 <Vector2Field U at node; step=4 ('UNLOADING2'); frame=0 >,
 <Tensor4Field S at node; step=4 ('UNLOADING2'); frame=1 >,
 <Vector2Field U at node; step=4 ('UNLOADING2'); frame=1 >]