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>