Orientation density functions


In [10]:
%matplotlib inline

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from simmit import smartplus as sim
from simmit import identify as iden
import os

dir = os.path.dirname(os.path.realpath('__file__'))

In this Python Notebook we will show how to properly run a simulation of a composite material, providing the ODF (orientation density function) of the reinforcments.

Such identification procedure require:

  1. Proper ODF peak data
  2. Proper composite properties
  3. A proper numerical model (here a composite model for laminate constitutive model)

In [12]:
x = np.arange(0,182,2)
path_data = dir + '/data/'
peak_file = 'Npeaks0.dat'

y = sim.get_densities(x, path_data, peak_file, False)
fig = plt.figure()

plt.grid(True)
plt.plot(x,y, c='black')


Out[12]:
[<matplotlib.lines.Line2D at 0x113b66590>]

In the previous graph we can see a multi-peak ODF (peaks are modeled using PEARSONVII functions). It actually represent quite well the microstructure of injected plates. The next step is to discretize the ODF into phases. The file containing the initial 2-phase microstructure contains the following informations


In [3]:
NPhases_file = dir + '/data/Nellipsoids0.dat'
NPhases = pd.read_csv(NPhases_file, delimiter=r'\s+', index_col=False, engine='python')
NPhases[::]


Out[3]:
Number Coatingof umat save c psi_mat theta_mat phi_mat a1 a2 a3 psi_geom theta_geom phi_geom nprops nstatev props
0 1 0 ELISO 1 0.8 0.0 0.0 0.0 1 1 1 0.0 0.0 0.0 3 1 3000
1 0 0 ELISO 1 0.2 0.0 0.0 0.0 50 1 1 0.0 0.0 0.0 3 1 70000

In [4]:
#Number_of_parameters
n_param = 6
#Number_of_consts
n_consts = 0
#Number_of_files
nfiles = 1

#Number_of_generations
ngen =  200
#Aleatory/Mesh space population : 0=mesh 1=meshlimit 2=random 3=defined
aleaspace = 2
#Space or aleatory population : apop in case of aleatory, spop in case of mesh
apop = 200
#Number of "doped" individual
ngboys = 1
#Max population per subgeneration
maxpop = 50
#stationnarity condition
station_nb = 20

path_data = dir + '/data'
path_keys = dir + '/keys'
path_results = dir + '/results'
outputfile = 'id_params.txt'
materialfile = 'material.dat'
simul_type = 'ODF'

iden.identification(simul_type,n_param,n_consts,nfiles,ngen,aleaspace,apop,ngboys,maxpop,station_nb,path_data,path_keys,path_results,materialfile,outputfile)


Out[4]:
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 0 ELISO 1 0.800000 0 0 0 1 1 1 0 0 0 3 1 3000
1 1 0 ELISO 1 0.024981 0 0 0 50 1 1 0 0 0 3 1 70000
2 2 0 ELISO 1 0.021216 5 0 0 50 1 1 5 0 0 3 1 70000
3 3 0 ELISO 1 0.014585 10 0 0 50 1 1 10 0 0 3 1 70000
4 4 0 ELISO 1 0.009644 15 0 0 50 1 1 15 0 0 3 1 70000

In [5]:
#Plot the concentration and the angle
c, angle = np.loadtxt(NPhases_file, usecols=(4,5), skiprows=2, unpack=True)

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()

# the histogram of the data
xs = np.arange(0,180,5)
rects1 = ax1.bar(xs, c, width=5, color='r', align='center')

ax2.plot(x, y, 'b-')

ax1.set_xlabel('X data')
ax1.set_ylabel('Y1 data', color='g')
ax2.set_ylabel('Y2 data', color='b')

ax1.set_ylim([0,0.025])
ax2.set_ylim([0,0.25])

plt.show()

#plt.grid(True)
#plt.plot(angle,c, c='black')
plt.show()



In [6]:
#Run the simulation
pathfile = 'path.txt'
nphases = 37 #The number of phases
num_file = 1 #The num of the file that contains the subphases
props = np.array([nphases, num_file, int1, int2])
outputfile = 'results_MTN.txt'
sim.solver(umat_name, props, nstatev, psi_rve, theta_rve, phi_rve, path_data, path_results, pathfile, outputfile)

fig = plt.figure()
outputfile_macro =  dir + '/' + path_results + '/results_MTN_global-0.txt'
e11, e22, e33, e12, e13, e23, s11, s22, s33, s12, s13, s23 = np.loadtxt(outputfile_macro, usecols=(8,9,10,11,12,13,14,15,16,17,18,19), unpack=True)
plt.grid(True)
plt.plot(e11,s11, c='black')

for i in range(8,12):
    outputfile_micro =  dir + '/' + path_results + '/results_MTN_global-0-' + str(i) + '.txt'
    e11, e22, e33, e12, e13, e23, s11, s22, s33, s12, s13, s23 = np.loadtxt(outputfile_micro, usecols=(8,9,10,11,12,13,14,15,16,17,18,19), unpack=True)
    plt.grid(True)
    plt.plot(e11,s11, c='red')

plt.xlabel('Strain')
plt.ylabel('Stress (MPa)')


plt.show()


---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-6-c7c8c1578bd9> in <module>()
      5 props = np.array([nphases, num_file, int1, int2])
      6 outputfile = 'results_MTN.txt'
----> 7 sim.solver(umat_name, props, nstatev, psi_rve, theta_rve, phi_rve, path_data, path_results, pathfile, outputfile)
      8 
      9 fig = plt.figure()

RuntimeError: Mat::operator(): index out of bounds

In [ ]: