CompMod2: Compartmentalized models with Argiope

Table of contents

Notebooks

Table of contents:

Moduli

Files in this folder:

1. local_settings.py
2. local_settings.py.example
3. model.pdf

Table of contents:

Note

This notebook can be downloaded here: DP_steel.ipynb

Dependency of the Young’s modulus to plastic strain in DP steels: a consequence of heterogeneity ?

Authors : Ludovic Charleux, Laurent Tabourot, Emile Roux, Moustapha Issack Farah, Laurent Bizet

DP980 steel tensile test simulation using a compartmentalized model.

Packages and various functions

%load_ext autoreload
%autoreload 2
%matplotlib nbagg

import compmod2 as cp2
import argiope as ag
import pandas as pd
import numpy as np
import inspect, os, local_settings, scipy
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d as mpl3d
import matplotlib.colors as colors
import matplotlib as mpl

# USEFUL FUNCTIONS
def create_dir(path):
    try:
        os.mkdir(path)
    except:
        pass
# SETTINGS
workdir   = "_workdir/"
outputdir = "outputs/"

label   = "RVE"
create_dir(workdir)
create_dir(workdir + outputdir)
create_dir(workdir + "reports/")

Model creation

The mode relies on two dedicated packages:

• CompMod2: compartmentalized model dedicated package.
• Argiope: FEM data/model management package.

Each parameter of the simulation can be tuned in this section and run into Abaqus.

shape = np.array([10, 10, 10]) # Number of elements in each direction for the RVE
Ne = shape.prod()

def element_map(mesh):
    """
    Element type definition.
    """
    mesh.elements.loc[:, ("type", "solver", "") ] = "C3D8R"
    return mesh

def material_map(mesh):
    """
    Material name definition
    """
    mesh.elements.materials = ["mat{0}".format(i) for i in mesh.elements.index]
    return mesh

def make_dist(k1, k2, l1, l2, w1):
    """
    Double Weibull distribution.
    """
    w2 = 1. - w1
    d1 = cp2.distributions.Weibull(k = k1, l = l1)
    d2 = cp2.distributions.Weibull(k = k2, l = l2)
    dist = cp2.distributions.CompositeDistribution(dists = [d1, d2], weights = [w1, w2])
    return dist


dist = make_dist(k1 = 1.64,    k2 = 3.17,
                 l1 = 4.16e-3, l2 = 3.42e-2,
                 w1 = 8.67e-1)

# Distribution Discretization Procedure (DDP)
xt, x = dist.discretize(Ne, xmax = 1.)

# random compartment element mapping (CEM)
RVE = np.arange(Ne)
np.random.shuffle(RVE) # Random RVE
x = x[RVE]


# Generation of the material of each compartment
E = 213.e3 # Young's modulus (MPa)
nu = .3    # Poisson ratio
materials = [ag.materials.ElasticPerfectlyPlastic(
                                 label = "mat{0}".format(i+1),
                                 young_modulus = E,
                                 poisson_ratio = nu,
                                 yield_stress = x[i] * E) # Spectification of the yield_stress for each compartment
                                 for i in range(Ne)]

# Definition of the laoding path aplly to the RVE
disp_values = np.linspace(0., .1, 10)[1:] # 9 Loading-unloading-reloading cycles
steps = []
for i in range(len(disp_values)):
    steps += [cp2.models.RVEStep(name = "loading{0}".format(i),
                                cx = ("disp", disp_values[i]),
                                field_output_frequency = 100),
             cp2.models.RVEStep(name = "unloading{0}".format(i),
                                cx = ("force", 0.),
                                field_output_frequency = 100),
             cp2.models.RVEStep(name = "reloading{0}".format(i),
                                cx = ("disp", disp_values[i]),
                                field_output_frequency = 100),]

sample = cp2.models.RVESample(shape = shape,
                              element_map = element_map,
                              material_map = material_map)

model  = cp2.models.RVETest(label = label,
                           parts = {"sample":sample},
                           steps = steps,
                           materials = materials,
                           solver = "abaqus",
                           solver_path = local_settings.ABAQUS_PATH,
                           workdir = workdir,
                           verbose = True)

Simulation Execusion

if False: # If True, the simulation will be run in Abaqus. Check load_settings.py before running this cell.
    model.write_input()
    model.run_simulation()
    model.postproc()
    model.save(workdir + label + ".pckl")

Simulation post-processing

All output data

model = ag.utils.load(workdir + label + ".pckl")
data = model.data["history"]
data.head()

	areas			dimensions			disp			energies	...	strains					stress			time	volume
	A1	A2	A3	L1	L2	L3	U1	U2	U3	We	...	E22	E33	LE11	LE22	LE33	S11	S22	S33	t	V
frame
0	0.999991	0.999991	0.999991	1.000000	1.000000	1.000000	0.000000	0.000000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.0	0.0	0.00	0.999991
1	0.999935	1.000080	1.000080	1.000111	0.999967	0.999967	0.000111	-0.000033	-0.000033	0.001312	...	-0.000033	-0.000033	0.000111	-0.000033	-0.000033	23.638306	0.0	0.0	0.01	1.000046
2	0.999873	1.000162	1.000162	1.000222	0.999933	0.999933	0.000222	-0.000067	-0.000067	0.005231	...	-0.000067	-0.000067	0.000222	-0.000067	-0.000067	47.157815	0.0	0.0	0.02	1.000095
3	0.999793	1.000226	1.000226	1.000333	0.999900	0.999900	0.000333	-0.000100	-0.000100	0.011713	...	-0.000100	-0.000100	0.000333	-0.000100	-0.000100	70.501260	0.0	0.0	0.03	1.000126
4	0.999727	1.000305	1.000305	1.000444	0.999866	0.999866	0.000444	-0.000134	-0.000134	0.020695	...	-0.000134	-0.000134	0.000444	-0.000134	-0.000134	93.606420	0.0	0.0	0.04	1.000171

5 rows × 33 columns

True strain vs. true stress curve

plt.figure()
for i in data.step.s.unique():
    step = data[data.step.s == i]
    color = "rgb"[i%3]
    plt.plot(step.strains.LE11, step.stress.S11, "-" + color,
             label = "step {0}".format(i))
plt.xlabel(r"Tensile True Strain, $\varepsilon$")
plt.ylabel(r"Tensile True Stress, $\sigma$")
plt.grid()
plt.show()

<IPython.core.display.Javascript object>

Moduli calculation

base on the protocol described in : Chen, Z., Gandhi, U., Lee, J., Wagoner, R.H., 2016b. Variation and consistency of young’s 15 modulus in steel. J. Mater. Process. Technol. 227, 227–243. (https://doi.org/10.1016/j.jmatprotec.2015.08.024)

moduli

	epsmax	epsp	mic	E1wag	E2wag	E3wag	E4wag	Ecwag
7	0.011050	0.006428	1.0	209090.338143	157098.837833	209313.102413	156177.948216	179912.149766
8	0.021979	0.016629	1.0	208178.024244	144962.584673	208781.412172	143510.235925	171676.021204
9	0.032790	0.027000	1.0	207649.990947	137960.719408	208267.703926	136193.294097	166814.630150
10	0.043485	0.037331	1.0	207266.591647	132680.462692	208033.880672	130698.949517	162958.952609
11	0.054067	0.047586	1.0	206973.833463	128381.920661	207970.526822	126304.004381	159688.966127
12	0.064539	0.057776	1.0	206803.623750	125055.502364	207622.096375	123169.620041	157097.918161
13	0.074901	0.067881	1.0	206654.886057	122104.666257	207302.743406	119703.357614	154817.595903
14	0.085158	0.077893	1.0	206524.706443	119377.042277	207092.278577	117499.145023	152742.003026
15	0.095310	0.087811	1.0	206409.660893	116858.180686	206808.821996	114622.127271	150811.872560

def modulus(strain, stress, low, high, order = 1):
    smax = stress.max()
    loc = (stress <= high * smax) & (stress >= low * smax)
    stress = stress[loc]
    strain = strain[loc]
    E = np.poly1d(np.polyfit(strain, stress, order)).deriv(1)(smax)
    return E

steps = data.step.s.unique()
moduli = pd.DataFrame(columns = ["epsmax", "epsp", "mic",
                              "E1wag", "E2wag", "E3wag", "E4wag", ])
trash, nc = moduli.shape
Ncycle = int(len(steps)/3)
m=1
for i in range(Ncycle):
    U = data[data.step.s == 3*i+1]
    R = data[data.step.s == 3*i+2]
    Eu = U.strains.LE11.values
    Su = U.stress.S11.values
    Er = R.strains.LE11.values
    Sr = R.stress.S11.values
    moduli.loc[m*nc+i, "mic"] = m
    moduli.loc[m*nc+i, "epsmax"] = Eu.max()
    moduli.loc[m*nc+i, "epsp"]  = Eu.min()
    moduli.loc[m*nc+i, "E1wag"] = modulus(Eu, Su, .7, 1.)
    moduli.loc[m*nc+i, "E2wag"] = modulus(Eu, Su, 0., .4)
    moduli.loc[m*nc+i, "E3wag"] = modulus(Er, Sr, 0., .3)
    moduli.loc[m*nc+i, "E4wag"] = modulus(Er, Sr, .6, 1.)
    moduli.loc[m*nc+i, "Ecwag"] = Su.max() / (Eu.max() - Eu.min())

moduli = moduli.applymap(float)

markers = "osv<*"
colors = "rgbmy"
mod = ["E1wag", "E2wag", "E3wag", "E4wag", "Ecwag"]
fig = plt.figure()
for i in range(len(mod)):
    plt.plot(moduli.epsmax * 100., moduli[mod[i]].values / 1000., markers[i] + "-", label = mod[i])
plt.grid()
plt.title("Moduli: Wagoner protocol")
plt.legend(loc = "best")
plt.xlabel("Log. (True) Strain, $LE_{11}$ [%]")
plt.ylabel("Moduli, $E$ [GPa]")
plt.show()

<IPython.core.display.Javascript object>

Deformed RVE state : 3D plot

step_label = "loading8"
frame = -1
field_label = "LE"
component_label = "v11"

# FIELDS MANAGEMENT
fdata = model.parts["sample"].mesh.fields_metadata()
fields = model.parts["sample"].mesh.fields
F_id = fdata[(fdata.step_label == step_label) & (fdata.label == field_label)].sort_values("frame").index[frame]
F = fields[F_id].data[component_label]
fdata.head()

	frame	frame_value	label	part	position	step_label	step_num
0	0	0	EE	ISAMPLE	element	loading0	0
1	0	0	LE	ISAMPLE	element	loading0	0
2	0	0	PE	ISAMPLE	element	loading0	0
3	0	0	S	ISAMPLE	element	loading0	0
4	0	0	U	ISAMPLE	node	loading0	0

# PATCHES
vertices, emap = model.get_Poly3DCollection(deformed = True,
                                            step_label = step_label,
                                            frame = frame,
                                            displacement_factor = 1.)

collection = mpl3d.art3d.Poly3DCollection(vertices)
collection.set_array(F.loc[emap] * 100)
collection.set_linewidth(.5)
collection.set_edgecolor("black")
collection.set_cmap(mpl.cm.jet)
collection.set_clim(0., 20.)

# FIGURE
fig = plt.figure()
ax1 = fig.add_subplot(1, 1, 1, projection='3d', aspect = "equal")
ax1.add_collection3d(collection)
#ax1.axis("off")
ax1.set_xlabel(r'loading direction, $\vec x$ axis')
ax1.set_ylabel(r'$\vec y$ axis')
ax1.set_zlabel(r'$\vec z$ axis')
offset = .05
ax1.set_xlim3d(-offset, 1. + offset)
ax1.set_ylim3d(-offset, 1. + offset)
ax1.set_zlim3d(-offset, 1. + offset)
ax1.set_xticks([0, 1])
ax1.set_yticks([0, 1])
ax1.set_zticks([0, 1])
cbar = plt.colorbar(collection)
cbar.set_label(r"True Strain Tensor, $LE_{11}$ [$\%$]")
cbar.set_ticks(np.linspace(0., 20, 6))
plt.savefig("model.pdf")
plt.show()

<IPython.core.display.Javascript object>

Rve

Files in this folder:

1. Makefile
2. conf.py
3. index.rst
4. local_settings.py
5. local_settings.py.example
6. model.pdf
7. notebooks_index.template

Table of contents:

Note

This notebook can be downloaded here: RVETest.ipynb

RVE Test

%load_ext autoreload
%autoreload 2

import compmod2 as cp2
import argiope as ag
import pandas as pd
import numpy as np
from string import Template
import inspect, os, local_settings, scipy
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d as mpl3d
import matplotlib.colors as colors
import matplotlib as mpl
%matplotlib nbagg

# USEFUL FUNCTIONS
def create_dir(path):
    try:
        os.mkdir(path)
    except:
        pass
# SETTINGS
workdir   = "_workdir/"
outputdir = "outputs/"
label   = "RVE"
create_dir(workdir)
create_dir(workdir + outputdir)

Model settings

%%time
shape = np.array([4, 4, 4])
Ne = shape.prod()

def element_map(mesh):
    mesh.elements.loc[:, ("type", "solver", "") ] = "C3D8R"
    return mesh


def material_map(mesh):
    mesh.elements.materials = ["mat{0}".format(i) for i in mesh.elements.index]
    return mesh

dist = cp2.distributions.Weibull(k = 1., l = 1.e-2)
xt, x = dist.discretize(Ne, xmax = 1.)
np.random.shuffle(x)
E = 210.e3
nu = .3
materials = [ag.materials.ElasticPerfectlyPlastic(
                                 label = "mat{0}".format(i+1),
                                 young_modulus = E,
                                 poisson_ratio = nu,
                                 yield_stress = x[i] * E)
                                 for i in range(Ne)]

steps = [cp2.models.RVEStep(name = "loading1",
                            cx = ("disp", .1),
                            field_output_frequency = 100),
         cp2.models.RVEStep(name = "unloading1",
                            cx = ("force", 0.),
                            field_output_frequency = 100),
         cp2.models.RVEStep(name = "loading2",
                            cx = ("disp", .1),
                            field_output_frequency = 100),]

sample = cp2.models.RVESample(shape = shape,
                              element_map = element_map,
                              material_map = material_map)
model = cp2.models.RVETest(label = label,
                           parts = {"sample":sample},
                           steps = steps,
                           materials = materials,
                           solver = "abaqus",
                           solver_path = local_settings.ABAQUS_PATH,
                           workdir = workdir,
                           verbose = True)

CPU times: user 10.5 ms, sys: 0 ns, total: 10.5 ms
Wall time: 10.5 ms

Solving

if True:
    model.write_input()
    #model.run_simulation()
    model.postproc()
    model.save(workdir + label + ".pckl")

#### POST-PROCESSING "RVE" USING POST-PROCESSOR "ABAQUS"
     Abaqus License Manager checked out the following license(s):
     "cae" release 6.13 from flex2-symme.univ-savoie.fr
     <7 out of 9 licenses remain available>.
  => POST-PROCESSED RVE: DURATION = 1.68s >

Post-Processing

model = ag.utils.load(workdir + label + ".pckl")
data = model.data["history"]
data.head()

	areas			dimensions			disp			energies	...	strains					stress			time	volume
	A1	A2	A3	L1	L2	L3	U1	U2	U3	We	...	E22	E33	LE11	LE22	LE33	S11	S22	S33	t	V
frame
0	1.000000	1.000000	1.000000	1.000	1.000000	1.000000	0.000	0.000000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	NaN	NaN	NaN	0.00	1.000000
1	0.999375	1.000687	1.000687	1.001	0.999688	0.999687	0.001	-0.000312	-0.000313	0.095873	...	-0.000312	-0.000313	0.001000	-0.000312	-0.000313	NaN	NaN	NaN	0.01	1.000374
2	0.998707	1.001349	1.001354	1.002	0.999356	0.999351	0.002	-0.000644	-0.000649	0.353291	...	-0.000644	-0.000649	0.001998	-0.000645	-0.000650	NaN	NaN	NaN	0.02	1.000704
3	0.998000	1.001992	1.002001	1.003	0.999004	0.998995	0.003	-0.000996	-0.001005	0.732349	...	-0.000996	-0.001005	0.002996	-0.000996	-0.001005	NaN	NaN	NaN	0.03	1.000994
4	0.997261	1.002617	1.002631	1.004	0.998636	0.998623	0.004	-0.001364	-0.001377	1.199574	...	-0.001364	-0.001377	0.003992	-0.001365	-0.001378	NaN	NaN	NaN	0.04	1.001250

5 rows × 33 columns

plt.figure()
for i in data.step.s.unique():
    step = data[data.step.s == i]
    plt.plot(step.strains.LE11, step.stress.S11, label = "step {0}".format(i))
plt.xlabel("Log. strain, $E$")
plt.ylabel("Cauchy Stress, $\sigma$")
plt.grid()
plt.legend(loc = "best")
plt.show()

<IPython.core.display.Javascript object>

plt.figure()
for i in data.step.s.unique():
    step = data[data.step.s == i]
    plt.plot(step.strains.LE11, step.energies.Wp, label = "step {0}".format(i))
plt.xlabel("Log. strain, $E$")
plt.ylabel("Plastic Dissipated Energy, $W_p$")
plt.grid()
plt.legend(loc = "best")
plt.show()

RVE plot

step_label = "loading1"
frame = -1
field_label = "LE"
component_label = "v11"

# FIELDS MANAGEMENT
fdata = model.parts["sample"].mesh.fields_metadata()
fields = model.parts["sample"].mesh.fields
F_id = fdata[(fdata.step_label == step_label) & (fdata.label == field_label)].sort_values("frame").index[frame]
F = fields[F_id].data[component_label]
fdata

	frame	frame_value	label	part	position	step_label	step_num
0	0	0	EE	ISAMPLE	element	loading1	0
1	0	0	LE	ISAMPLE	element	loading1	0
2	0	0	PE	ISAMPLE	element	loading1	0
3	0	0	S	ISAMPLE	element	loading1	0
4	0	0	U	ISAMPLE	node	loading1	0
5	1	1	EE	ISAMPLE	element	loading1	0
6	1	1	LE	ISAMPLE	element	loading1	0
7	1	1	PE	ISAMPLE	element	loading1	0
8	1	1	S	ISAMPLE	element	loading1	0
9	1	1	U	ISAMPLE	node	loading1	0
10	0	0	EE	ISAMPLE	element	unloading1	1
11	0	0	LE	ISAMPLE	element	unloading1	1
12	0	0	PE	ISAMPLE	element	unloading1	1
13	0	0	S	ISAMPLE	element	unloading1	1
14	0	0	U	ISAMPLE	node	unloading1	1
15	1	1	EE	ISAMPLE	element	unloading1	1
16	1	1	LE	ISAMPLE	element	unloading1	1
17	1	1	PE	ISAMPLE	element	unloading1	1
18	1	1	S	ISAMPLE	element	unloading1	1
19	1	1	U	ISAMPLE	node	unloading1	1
20	0	0	EE	ISAMPLE	element	loading2	2
21	0	0	LE	ISAMPLE	element	loading2	2
22	0	0	PE	ISAMPLE	element	loading2	2
23	0	0	S	ISAMPLE	element	loading2	2
24	0	0	U	ISAMPLE	node	loading2	2
25	1	1	EE	ISAMPLE	element	loading2	2
26	1	1	LE	ISAMPLE	element	loading2	2
27	1	1	PE	ISAMPLE	element	loading2	2
28	1	1	S	ISAMPLE	element	loading2	2
29	1	1	U	ISAMPLE	node	loading2	2

# PATCHES
vertices, emap = model.get_Poly3DCollection(deformed = True,
                                            step_label = step_label,
                                            frame = frame,
                                            displacement_factor = 1.)

collection = mpl3d.art3d.Poly3DCollection(vertices)
collection.set_array(F.loc[emap] * 100)
collection.set_linewidth(.5)
collection.set_edgecolor("black")
collection.set_cmap(mpl.cm.jet)
collection.set_clim(0., 20.)

# FIGURE
fig = plt.figure()
ax1 = fig.add_subplot(1, 1, 1, projection='3d', aspect = "equal")
ax1.add_collection3d(collection)
#ax1.axis("off")
ax1.set_xlabel(r'loading direction, $\vec x$ axis')
ax1.set_ylabel(r'$\vec y$ axis')
ax1.set_zlabel(r'$\vec z$ axis')
offset = .05
ax1.set_xlim3d(-offset, 1. + offset)
ax1.set_ylim3d(-offset, 1. + offset)
ax1.set_zlim3d(-offset, 1. + offset)
ax1.set_xticks([0, 1])
ax1.set_yticks([0, 1])
ax1.set_zticks([0, 1])
cbar = plt.colorbar(collection)
cbar.set_label(r"True Strain Tensor, $LE_{11}$ [$\%$]")
cbar.set_ticks(np.linspace(0., 20, 6))
plt.savefig("model.pdf")
plt.show()

<IPython.core.display.Javascript object>

Periodic BC checking

U_id = fdata[(fdata.step_label == step_label) & (fdata.label == "U")].sort_values("frame").index[frame]
U = fields[U_id].data
nodes = model.parts["sample"].mesh.nodes

def check_direction(direction = "x"):
    """
    Calculates the difference of the displacements on 2 opposite faces that are normal to a given direction.
    """
    F0 = nodes[(nodes.sets.control == False) & (nodes.coords[direction] == 0.)].coords.sort_values(["x", "y", "z"])
    F1 = nodes[(nodes.sets.control == False) & (nodes.coords[direction] == 1.)].coords.sort_values(["x", "y", "z"])
    U0 = U.loc[F0.index].reset_index()
    U1 = U.loc[F1.index].reset_index()
    return U1 - U0

check_direction("x").describe()

	node	v1	v2	v3
count	25.0	2.500000e+01	25.0	25.0
mean	4.0	1.000000e-01	0.0	0.0
std	0.0	2.665862e-09	0.0	0.0
min	4.0	1.000000e-01	0.0	0.0
25%	4.0	1.000000e-01	0.0	0.0
50%	4.0	1.000000e-01	0.0	0.0
75%	4.0	1.000000e-01	0.0	0.0
max	4.0	1.000000e-01	0.0	0.0

check_direction("y").describe()

	node	v1	v2	v3
count	25.0	25.0	2.500000e+01	25.0
mean	20.0	0.0	-4.163784e-02	0.0
std	0.0	0.0	7.457459e-10	0.0
min	20.0	0.0	-4.163784e-02	0.0
25%	20.0	0.0	-4.163784e-02	0.0
50%	20.0	0.0	-4.163784e-02	0.0
75%	20.0	0.0	-4.163784e-02	0.0
max	20.0	0.0	-4.163784e-02	0.0

check_direction("z").describe()

	node	v1	v2	v3
count	25.0	25.0	25.0	2.500000e+01
mean	100.0	0.0	0.0	-4.791145e-02
std	0.0	0.0	0.0	1.103872e-09
min	100.0	0.0	0.0	-4.791145e-02
25%	100.0	0.0	0.0	-4.791145e-02
50%	100.0	0.0	0.0	-4.791145e-02
75%	100.0	0.0	0.0	-4.791145e-02
max	100.0	0.0	0.0	-4.791145e-02

The boundary conditions are perfectly matched.