Tutorials
Please note: Make sure you activate your environment before running any python code that invoke Simcoon.
The python script and input files for this tutoial can be found here:
conda activate scientific
To be able to follow these tutorials it is much better to install the ipython environment and a few more libraries:
conda install ipython
conda install ipykernel
python -m ipykernel install --user --name=scientific
conda install -c conda-forge matplolib matplotlib-inline pandas
Composite simulation with mean field methods
Please note: Be careful with the data type while you pass numpy arrays. Make sure you define them as dtype=’float’ unless specified, since they will be interpreted as float by Simcoon, i.e.: x = np.array([1, 2, 3, 4, 5], dtype=’float’)
In this tutorial we will study the evolution of the mechanical properties of a composite material considering spherical reinforcement (sew figure below)
Elastic properties will be evaluated depending on the volume fraction of reinforcement (up to 50% volume fraction).
The following elastic properties for the matrix is considered: \(E = 2250\) MPa, \(\nu = 0.19\). The following elastic properties for the reinforcement is considered: \(E = 2250\) MPa, \(\nu = 0.19\)
A: Effective properties of the composite with 20% of reinforcement
The first thing we want to do is to add a file in a ‘data’ folder, named ‘Nellipsoids0.dat’. This file can be downloaded here:
Number | Coatingof | umat | save | c | psi_mat | theta_mat | phi_mat | a1 | a2 | a3 | psi_geom | theta_geom | phi_geom | nprops | nstatev | props | ||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | ELISO | 1 | 0.8 | 0. | 0. | 0. | 1 | 1 | 1 | 0. | 0. | 0. | 3 | 1 | 2250 | 0.19 | 0. |
1 | 0 | ELISO | 1 | 0.2 | 0. | 0. | 0. | 1 | 1 | 1 | 0. | 0. | 0. | 3 | 1 | 73000 | 0.19 | 0. |
import numpy as np
from simcoon import simmit as sim
Next we shall inform the number of internal state variables (if any) at the macroscopic level and the material properties:
nstatev = 0 #None here
nphases = 2 #The number of phases
num_file = 0 #The num of the file that contains the subphases
int1 = 50 #Number of integration points in the long axis
int2 = 50 #Number of integration points in the lat axis
n_matrix = 0 #phase number for the matrix
props = np.array([nphases, num_file, int1, int2, n_matrix], dtype='float')
There is a possibility to consider a misorientation between the test frame of reference and the composite material orientation, considering Euler angles with the z-x-z. The stiffness tensor will be returned in the test frame of reference. Also, we shall inform which micrmechanical scheme will be used. Enter umat_name = ‘MIMTN’ for a Mori-Tanaka scheme and umat_name = ‘MISCN’ for a self-consistent scheme
psi_rve = 0.
theta_rve = 0.
phi_rve = 0.
umat_name = 'MIMTN'
We are now ready to run the simulation, get the effective isotropic elastic properties (since both phases are isotropic and the reinforcements are spherical) and print them:
L = sim.L_eff(umat_name, props, nstatev, psi_rve, theta_rve, phi_rve)
p = sim.L_iso_props(L)
print(p)
The result is a python numpy array containing the material parameters (here two, the Young’s modulus and the Poisson ratio):
[3.29316019e+03 1.91800053e-01]
Notice: Note that despite the Poisson ratio is the same between the two materials, the stiffness mismatch between the two phases lead to a different effective Poisson ratio.
B: Effective properties of the composite as a function of the reinforcements
02: Cubic symmetry and directional stiffness
In this tutorial we will plot the directional stiffness of a cubic material. Make sure first that you have matplotlib installed:
conda install -c conda-forge matplolib
We shall first import the libraries required, i.e. lumpy, matplotlib and simcoon. We also indicate to use \(\LaTeX\) into the matplotlib plots:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm, colors
from simcoon import simmit as sim
We explore next the possible direction in the three-dimensional space using a unit vector \(\vec{n}\), and the result is stored into a numpy array:
\(\vec{n} = \left( \begin{array}{c} \sin (\theta) \cos (\phi) \\ \sin (\theta) \sin (\phi) \\ \cos (\theta) \end{array} \right)\)
phi = np.linspace(0,2*np.pi, 128) # the angle of the projection in the xy-plane
theta = np.linspace(0, np.pi, 128).reshape(128,1) # the angle from the polar axis, ie the polar angle
n_1 = np.sin(theta)*np.cos(phi)
n_2 = np.sin(theta)*np.sin(phi)
n_3 = np.cos(theta)*np.ones(128)
n = np.array([n_1*n_1, n_2*n_2, n_3*n_3, n_1*n_2, n_1*n_3, n_2*n_3]).transpose(1,2,0).reshape(128,128,1,6)
Next we are using Simcoon to determine the elastic stiffness matrix \(\mathcal{L}\) from the cubic stiffness components \((C_{11}, C_{12}, C_{44})\) and to invert it to obtain the compliance matrix \(\mathcal{M}\):
C11 = 185000.
C12 = 158000.
C44 = 39700.
L = sim.L_cubic(C11, C12, C44, 'Cii')
M = np.linalg.inv(L)
We write now a numpy array that performs the operation \(\vec{n} \cdot \mathcal{M} \cdot \vec{n}\) and invert it to obtain the elastic stiffness for all the directions of the vector \(\vec{n}\). The last operation is to get the components \((x,y,z)\) of this stiffness to be able to plot a surface plot:
S = (n@M@n.reshape(128,128,6,1)).reshape(128,128)
E = (1./S)
x = E*n_1
y = E*n_2
z = E*n_3
The last part consist in plotting the results into a surface plot and get the nice following picture:
fig = plt.figure(figsize=plt.figaspect(1)) # Square figure
ax = fig.add_subplot(111, projection='3d')
#ax.plot_surface(x, y, z, cmap='hot',c=E)
norm = colors.Normalize(vmin = np.min(E), vmax = np.max(E), clip = False)
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, norm=norm, facecolors=cm.bone(norm(E)),linewidth=0, antialiased=False, shade=False)
#ax.set_xlim(0,20000)
#ax.set_ylim(0,20000)
#ax.set_zlim(0,20000)
ax.set_xlabel(r'$E_x$')
ax.set_ylabel(r'$E_y$')
ax.set_zlabel(r'$E_z$')
scalarmap = cm.ScalarMappable(cmap=plt.cm.bone, norm=norm)
#m.set_array([])
cbar = plt.colorbar(scalarmap)
#cbar.ax.set_yticklabels(['0','1','2','>3'])
cbar.set_label(r'directional stiffness $E$', rotation=270)
plt.show()